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ABSTRACT 
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This  dissertation  presents  a degree  of  freedom  or  information 
content  analysis  of  images  and  imaging  systems  in  the  context  of 
digital  image  processing.  As  such  it  represents  an  attempt  to 
quantify  the  number  of  truly  independent  samples  one  gathers  with 
imaging  devices. 

In  quantifying  the  degrees  of  freedom  of  an  imaging  system  it  is 
necessary  to  develop  an  appropriate  model.  In  this  work  the  imaging 
system  is  modeled  as  a linear  system  through  the  cont inuous -discret e 
imaging  equation.  The  associated  gram  matrix  is  then  employed  as 
an  aid  in  defining  the  system  degrees  of  freedom.  The  gram  matrix 
eigenvalues  are  shown  to  be  related  to  those  of  the  associated 
continuous -continuous  model  and  can  be  used  to  predict  the  discre- 
tized system  performance.  These  ideas  are  then  applied  to  the 
tomographic  or  projection  imaging  system,  and  result  in  the  ability 
to  predict  the  performance  of  this  system  by  indicating  where 
redundant  data  is  achieved,  and  the  best  ways  of  increasing  the 
degrees  of  freedom  with  a minimum  sample  increase. 

The  degrees  of  freedom  of  a sampled  image  itself  are  developed 
as  an  approximation  problem.  Here  bicubic  splines  with  variable 
knots  are  employed  in  an  attempt  to  answer  the  question  as  to  what 
extent  images  are  finitely  representable  in  the  context  of  a digital 

ii 


computer. 

Relatively  simple  algorithms  for  good  knot  placement  are  given, 
and  result  in  spline  approximations  that  achieve  significant  parameter 
reductions  at  acceptable  error  levels.  The  knots  themselves  are 
shown  to  be  useful  as  an  indicator  of  image  activity,  and  have 
potential  as  an  image  segmentation  device. 
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Chapter  1 


INTRODUCTION 

This  dissertation  is  concerned  with  the  concept  of  degrees  of 
freedom  (DOF)  or  information  content  of  images  and  imaging  systems 
arising  in  digital  image  processing.  This  concept  is  important,  as  it 
is  fundamental  to  problems  such  as  image  coding  and  image  restora- 
tion. In  coding  problems,  one  is  interested  in  the  transmission  of 
that  information  relevant  to  the  users'  needs  and  in  the  elimination  of 
irrelevant  data;  while  in  the  restoration  problems,  the  object  is  to  be 
restored  from  samples  of  a corrupted  image. 

These  two  examples  are  given  to  illustrate  that  the  subject  could 
be  approached  from  two  points  of  view.  Namely,  the  image  could  be 
treated  as  the  output  of  an  imaging  system  whose  characteristics  are 
known  and  thus  dependent  on  the  DOF  of  the  imaging  system;  or  the 
subject  could  be  considered  by  itself  as  in  the  coding  problem  where 
the  imaging  system  characteristics  are  either  secondary  in  import- 
ance or  considered  ideal. 

Fundamentally,  this  concept  of  degrees  of  freedom  can  be  viewed 
as  an  attempt  to  quantify  the  number  of  truly  independent  samples  of 
data  one  gathers  with  photographic  or  other  imaging  devices.  As 
image  sensor  technology  grows,  the  quantity  of  data  gathered 


increases,  and  it  becomes  reasonable  to  ask  what  the  true  increase  in 


| 

. 


i 

f, 


information  content  is  as  one  increases  image  samples.  This  is 
especially  important  in  medical  imaging  applications  where  an 
increase  in  the  quantity  of  data  gathered,  while  not  producing  a 
corresponding  increase  in  image  information,  subjects  the  patient  to 
an  unnecessary  increase  in  radiation  exposure.  Thus  in  designing 
imaging  systems  for  medical  applications  it  is  extremely  important 
that  the  information  content  of  the  imaging  system  be  quantized. 

Since  half  the  thrust  of  this  dissertation  will  be  towards  inform- 
ation content  in  imaging  systems  in  general,  and  to  the  tomographic, 
or  projection,  imaging  system  in  particular,  a mathematical  model 
for  imaging  will  be  necessary. 

1 . 1 Mathematical  Imaging  Models 

In  modelling  imaging  systems  an  assumption  often  made  is  that 
of  a linear  system.  While  it  is  not  true  that  every  imaging  system  is 
in  fact  linear,  this  assumption  is  useful  in  that  it  makes  the  analysis 
tractable  and  provides  reasonable  results.  Even  if  untrue,  the 
system  can  often  be  considered  linear  if  the  region  of  observation  is 
kept  small  enough.  Thus  imaging  systems  can  be  analyzed  by 
considering  them  in  terms  of  two-dimensional  linear  system  theory. 
In  applying  these  methods  to  imaging  systems  the  assumption  is 
made  that  an  image,  g,  is  related  to  the  original  object,  f,  by  a 
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superposition  integral  as  follows; 

g(x,y)  = JT  h(x,  y;5,  n)f(c,  T))dFdTl  . (1-1) 

R 

Here  h(x,  y;F,  T |)  represents  the  weighting  function  of  the  imaging 
system  and  R is  the  region  of  integration  over  the  input  coordinate 
system.  In  this  dissertation  the  assumptions  will  be  made  that 
h(x,  y;f,T1  ) is  continuous  in  x and  y,  bounded,  and  such  that 

[ITJ  |h(x,  y;P,  Tl)  |2  dxdydcdT)  <® 

Furthermore,  it  will  be  assumed  that  f(^  , T»  > will  also  be  bounded  and 

at  least  piecewise  continuous.  These  assumptions  will  imply  that 

g(x,  y)  will  be  continuous  in  x and  y.  This  assumption  on  the  part  of 

g(x,  y)  is  necessary  for  digital  processing  since  only  a sampled  , 

version  of  g will  be  dealt  with.  Thus  if  we  define  g.  to  be  g(x,  y) 

sampled  in  the  image  coordinate  plane  (x,  y)  and  similarly  for  h,  then 

the  vector  form  of  equation  (1-1),  the  continuous -discrete  imaging 

model,  is  obtained  [l-l,  1 — Z 1 

& = ff  h(p,  Tl)f(c,  D)d*dTl  + n (1-2) 

R 

where  ri  is  an  error  term  and  all  vectors  are  Nxl  columns.  Here  the 
aim,  as  in  any  reconstruction  process,  is  to  recover  f as  best  as 
possible,  in  some  sense,  knowing  h and  This  model  is  depicted 
in  Fig.  (1.  1). 

3 


Intuitively,  the  concept  of  degrees  of  freedom  (DOF)  fits  nicely 
into  this  model  as  can  be  seen  by  noting  that  f(Ff  tl)  js  defined  on  a 
continuum,  namely  R,  and  as  such  represents  a noncountably 
infinite  number  of  DOF.  The  image  g.  is  described  by  a finite 
number  of  samples  N,  and  at  best  represents  N DOF.  However,  the 
imaging  integral  equation  causes  reduction  of  this  number  due  to  the 
point  spread  function  (PSF)  blur  provided  by  h(5',  Tl).  Here  the  PSF 
is  described  by  an  N-vector  whose  elements  are  continuous  functions 
of  the  object  coordinate  system  (P,T).  One's  intuition  might  serve 
here  in  the  concept  that  the  greater  the  blur,  the  fewer  the  DOF. 
Thus,  if  h.(=-,  7) ) is  a "narrow"  function  the  DOF  would  be  greater 
than  if  h^r,  1)  ) were  a "broad"  function.  Because  h.C?,  T ) can  repre- 
sent  a space  variant  point  spread  function  (SVPSF),  we  can  further 
let  our  intuition  suggest  that  in  regions  in  which  the  object  is  in 
"better  focus"  (i.  e.,  narrower  PSF)  we  would  be  obtaining  greater 

DOF  than  in  regions  of  "poorer  focus"  (i.e.,  broader  or  greater  blur 
PSF). 

The  rationale  behind  the  separation  of  the  problem  into  two 
I 

classes  should  be  a little  clearer  now.  If  the  system  is  taken  to  be 
ideal,  then  £ is  simply  f(f,  7)  ) sampled  in  the  output  coordinate  plane 
and  we  are  confronted  with  the  problem  of  relating  the  DOF  of  a 
sampled  image  to  its  original  unsampled  version.  It  must  also  be 
considered  that  to  sensibly  discuss  sampled  images,  one  must  be 
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dealing  with  ones  that  are  band  limited.  If  f(§,  T1  ) is  not  band  limited, 
which  is  often  the  case,  we  can  assume  that  its  sampled  version  f_  is 
obtained  from  some  other  function  whose  band  limited  version 


coincides  with  f at  the  sample  points.  In  this  case,  h can  be  taken 
as  an  ideal  low  pass  filter. 

1 . 2 Research  Objectives 

The  objectives  of  this  research  will  be  to  develop  a measure  of 
the  DOF  of  imaging  systems  and  to  apply  this  to  two  areas  of  study. 

The  first  is  that  of  the  tomographic  or  projection  imaging  system. 

The  second  area  will  be  that  of  developing  "smart  sensors"  by  variably 
adapting  to  the  DOF  of  the  sensed  image.  With  respect  to  the  former 
projection  imaging  system  [l  -3]  , it  will  be  necessary  to  develop  a 
weighting  function  for  the  tomographic  imaging  system  and  its 
associated  continuous -discrete  model. 

By  quantifying  the  degrees  of  freedom  of  the  tomographic  imaging 
system  we  will,  for  the  first  time,  be  able  to  predict  the  resolution 
capability  of  the  system  for  large  numbers  of  samples.  It  will  also  be 
shown  that  the  resolution  limits  in  tomography  are  not  a function  of 
the  particular  reconstruction  algorithm  employed,  but  are  funda- 
mental to  the  process  itself. 

While  the  main  intent  in  the  first  area  of  research  is  to  investi- 
gate the  degrees  of  freedom  and  information  content  in  the  projection 
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imaging  process,  structure  will  be  shown  to  exist  that  makes  a 
linear  algebraic  solution  to  the  reconstruction  problem  feasible  for 
images  of  a rather  large  dimension.  This  linear  algebraic  algorithm 
will  be  developed,  and  the  results  presented  in  this  dissertation 
concerning  information  content  in  projection  imaging  will  be  obtained 
from  this  model. 

It  will  also  be  shown  that  the  difficulty  in  image  restoration 
arises  not  so  much  from  the  sampling  scheme  used,  but  is  innate  to 
the  original  continuous -continuous  model  of  linear  imaging  systems 
in  general.  This  will  be  accomplished  by  relating  the  eigenvalues  of 
the  gram  matrix  of  the  kernel  vector  h with  the  singular  values  of 
the  original  kernel  h(x,  y;?,  1]).  The  gram  matrix  will  be  discussed  in 
Ch.  3.  The  situation  where  h(x,  y ; ~ , T> ) is  unknown  will  also  be  treated. 
Here  the  problem  will  be  considered  as  a two  dimensional  approxi- 
mation problem  and  the  concept  of  an  "epsilon  degrees  of  freedom" 
will  be  developed.  By  this  it  is  meant  that  the  degrees  of  freedom 
of  an  image  at  a level  epsilon  will  be  the  minimum  number  of  func- 
tions needed  to  approximate  f(ff,  11  ) within  an  accuracy  of  epsilon 
assuming  a particular  metric. 

From  a "smart  sensor"  viewpoint  by  way  of  motivation,  if  we 
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consider  a sampled  image  consisting  of  N samples  that  could  be 
approximated  to  an  acceptable  error  by  a least  squares  polynomial  of 
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M variables  with  M , it  would  be  reasonable  to  say  that  this 

2 

image  had  less  than  N DOF  in  a least  squares  sense  using  a poly- 
nomial as  the  approximation  technique.  This  approach  is  taken  to 
circumvent  the  difficulty  of  associating  a finite  DOF  to  an  image 
defined  on  a continuum  which  obviously  has  an  uncountably  infinite 
DOF  if  we  desire  to  specify  that  image  exactly.  However,  if  we  are 
willing  to  accept  an  approximation  with  a small  but  nonzero  error, 
then  the  possibility  exists  in  quantifying  the  DOF  in  this  manner. 

In  effect  this  represents  an  attempt  to  "bridge  the  gap"  between 
the  continuous  domain  upon  which  images  are  defined  and  the  discrete 
point  sets  involved  in  digital  computations.  Shannon's  sampling 
theory  represents  one  method  for  reconstructing  uniformly  sampled 
band  limited  images,  but  by  taking  a more  general  approximation 
theoretic  approach  other  sampling  and  reconstruction  techniques  may 
be  used.  For  example,  the  desirability  of  nonuniform  or  adaptive 
sampling  can  be  illustrated  by  considering  an  image  that  contains 
high  frequency  information  in  a small  region  in  its  domain  of 
definition  with  the  remainder  containing  only  low  frequency  compon- 
ents. If  it  is  desired  to  reconstruct  this  sampled  image  using 
Shannon  sampling,  then  it  must  be  sampled  at  the  Nyquist  rate 
defined  by  that  small  high  frequency  zone.  Intuitively,  it  would  seem 
then  that  the  regions  of  low  frequency  content  are  over  sampled. 
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Realizing  that  we  can  never  collect  an  infinite  number  of  samples  for 
application  on  a digital  computer  and  that  some  error  always  results 
from  this  technique,  it  might  be  possible  to  achieve  a reasonable 
error  by  a finite  "Shannon  interpolation"  at  the  respective  Nyquist 
rates  in  each  of  those  regions.  More  formally,  in  one  dimension  if 
is  the  bandwidth  of  the  function  and  T is  the  length  of  the  interval 
of  observation  then  there  are  approximately  2 B T independent 
samples  from  this  function  in  the  interval  [o,  T]  . If  we  wish  to  say 
that  in  the  interval  [o,a],  a < T the  hignest  "effective"  frequency 
component  is  BT  and  in  [a,  T]  the  highest  "effective"  frequency 
component  is  Bj,  (B^  < B ),  then  the  number  of  independent  samples 
in  [0,Tl  becomes  2 [b Ta  + B^T-a)].  Since  Bj  < B it  follows  that 
^[B^a  + Bj  (T-a)l  < ZB^T.  Depending  on  the  ratio  of  B,^  to  Bj  it 
might  be  possible  to  reconstruct  the  function  to  reasonable  error  by 
this  approach  with  far  fewer  samples. 

To  illustrate  the  applicability  of  adaptive  processing  consider 
that  advances  in  charge  coupled  device  (CCD)  sensors  are  such  that 
some  preprocessing  within  the  sensor  itself  is  not  so  unrealistic, 
this  preprocessing  could  involve  some  evaluation  as  to  what  data 
constitutes  information  to  the  user  and  transmits  only  that  data 
relevant  to  the  users'  needs.  Surface  acoustic  wave  devices  (SAW), 
are  becoming  available  that  can  provide  a Fourier  transform  of  an 
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image  at  video  data  rates  so  that  it  is  possible  to  obtain  a sensor  that 
provides  a transform  of  the  image  as  its  output.  This  increase  in 
sensor  sophistication  coupled  with  the  ability  to  gather  large  quantities 
of  data,  the  ability  to  do  adaptive  sampling  or  some  other  more 
exotic  processing  to  get  at  the  real  information  content  in  the  data, 
may  provide  fruitful  results. 


Chapter  2 
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REVIEW  OF  THE  STATE  OF  THE  ARI 

2.  1 Tomographic  Literature  Review 

The  concept  of  degrees  of  freedom  (DOF)  has  arisen  in  the 
imaging  literature  in  an  attempt  to  quantify  the  number  of  truly 
independent  samples  of  data  one  gathers  with  photographic  or  other 
image  sensing  devices. 

In  obtaining  an  estimate,?,  for  f in  (1-2),  consideration  must  be 
taken  of  the  fact  that  the  PSF  kernel  vector,  h>  is  comprised  of 
kernel  functions  that  are  not  necessarily  independent.  Neglecting 
the  error  term  for  the  moment,  this  dependency  of  kernel  functions 
implies  that  some  of  the  output  samples  can  be  predicted  by  linear 
combinations  of  the  others,  are  thus  superfluous,  and  serve  to  reduce 
the  DOF  of  the  imaging  system.  The  error  term  aggravates  this 
situation  in  so  far  as,  if  an  output  element  can  be  predicted  by  a 
linear  combination  of  the  others  to  within  an  accuracy  better  than  the 
measurement  error,  that  measurement  adds  no  new  information  to 
what  is  already  known.  For  one  dimensional  systems  it  has  been 
suggested  that  the  Gramian  formed  from  the  PSF  N-vector  could 
provide  a quantitative  measure  of  the  DOF  available  in  an  imaging 
system.  The  entries  in  this  matrix  represent  the  correlation  or 
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overlap  of  each  PSF  with  its  neighbors.  If  we  consider  the  following 


one  dimensional  continuous -di  screte  model 


The  gram  matrix  [r]  is  given  by 

[r]  = f hr)h;|:(?)d* 

K 

where  all  vectors  are  considered  as  Nxl  column  vectors  and  * 
indicates  vector  conjugate  transpose.  The  rank  and/or  eigenvalue 
map  could  be  used  in  the  definition  of  the  DOF  [2-1,2-2,2-31.  In 
this  dissertation  we  pursue  this  development  in  two  dimensions  and 
apply  the  analysis  to  the  tomographic  or  projection  imaging  system 
as  an  imaging  system  example. 

Projection  imaging  or  three-dimensional  reconstruction  has  been 
the  subject  of  much  research  in  recent  years.  This  is  the  result  of 
important  applications  such  as  electron  microscopy  [2-4,2-51,  radio 
astronomy  [2-6,  2-7],  and  trans -axial  tomography  [2-8,  2-9,  2-10],  to 
name  a few.  There  has  been  a somewhat  intense  effort  in  this  last 
application  in  the  medical  community  with  the  proliferation  of 
articles,  reports,  conferences,  and  even  the  manufacture  of 
equipment.  While  the  mathematical  basis  upon  which  transaxial 
tomography  is  founded  is  quite  sound,  many  of  the  practioners 
of  the  technique  have  developed  a variety  of  image  reconstruction 
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approaches  which  are  seemingly  unrelated  and  controversial. 


Basically  the  technique  requires  computational  reconstruction  of 
an  image  and  researchers  have  developed  algorithms  which  can  be 
roughly  classified  as  follows:  convolutional  [2-6,2-71,  Fourier 
transform  [2-4]  , algebraic,  in  the  sense  of  ART  [2-91  , and  linear 
algebraic  as  described  by  Kashyap  [2-lo!  . The  convolutional  and 
Fourier  transform  algorithms  can  be  considered  as  closed-form 
algorithms  in  the  sense  that  a continuous  solution  is  first  obtained, 
and  then  discretized  for  implementation  on  a digital  computer. 

While  mathematically  elegant,  this  sequence  gives  little  insight  into 
the  errors  incurred  due  to  discretization  and  the  presence  of  noise 
or  measurement  errors  in  the  system.  Linear  algebraic  approaches 
somewhat  circumvent  this  situation,  but  have  not  received  as  much 
attention,  most  likely,  as  large  matrices  are  involved  for  images 
of  even  a moderate  dimension.  Nevertheless,  the  information 
content  in  a particular  sampling  geometry  can  be  determined  by 
considering  the  spectrum  of  the  point  spread  function  matrix 
describing  the  imaging  system.  This  is  not  to  say  that  l.near 
algebraic  methods  are  the  only  way  to  approach  the  information 
content  or  error  analysis  in  the  projection  imaging  system.  King 
and  Crowfher  [2-111  have  formulated  the  image  reconstruction 
problem  as  an  eigenvalue  problem  in  the  continuous  domain,  and 
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have  drawn  analogies  to  optical  data  processing  systems  in  analyzing 
the  effect  of  a finite  number  of  projections.  In  addition  to  this  effort. 

Smith,  Peters  and  Bates  [2-12]  have  considered  the  effect  of  a finite 
number  of  projections  in  the  reconstruction  process. 

2 . 2 Image  Approximation  Literature  Review 

The  concept  of  the  degrees  of  freedom  of  a sampled  image  as  an 
approximation  problem  arises -quite  naturally  in  the  context  of  image 
coding  by  transform  methods.  In  transform  image  coding  an  ortho- 
gonal transformation  is  performed  on  a sampled  image  matrix  and  a 
bandwidth  reduction  is  obtained  by  transmitting  only  those  transform 
coefficients  above  a certain  threshold  whose  level  is  consistent  with 
the  desired  error  [2-13],  The  large  bandwidth  reductions  reported 
are  due,  in  part,  to  the  compacting-of-image -energy  property  of 
the  orthogonal  transforms  employed.  However,  any  compacting  in 
the  transform  domain  is  at  the  expense  of  an  increased  dynamic  range 
in  the  transform  coefficients  because  of  the  conservation  of  energy 
inherent  to  all  orthogonal  transformations. 

Another  situation  is  the  application  of  the  singular  value 
decomposition  (SVD)  algorithm  [2-14]  to  the  sampled  image  matrix 
where  upon  the  number  of  degrees  of  freedom  can  be  equated  with  the 
number  of  effectively  nonzero  singular  values,  with  the  remaining 
parameters  describing  the  orthogonal  singular  vectors. 
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In  both  of  these  situations  we  are  dealing  with  a sampled  version 


of  the  image  only,  and  as  a result  the  degrees  of  freedom  could  be 
affected  by  the  sampling  method  used  and  thus  could  lead  to  mis- 
leading results.  The  point  to  be  made  here  is  that  the  degrees  of 
freedom  should  be  a characteristic  of  the  original  image  and 
reflected  in  the  sampled  image  only  by  our  inability  to  collect  an 
uncountably  infinite  number  of  samples  for  application  on  a computer. 

It  is  to  this  end  that  approximation  theory  is  directed.  In  fact 
the  Remes  algorithm  represents  a method  for  finding  the  best 
approximating  polynomial  (in  the  uniform  norm)  of  degree  Sntoa 
continuous  function  in  the  interval  [a,b]  in  terms  of  a sequence  of 
solutions  involving  finite  sample  point  sets  with  n+2  elements  each 
[2-1  5]  . Unfortunately,  no  equivalent  algorithm  is  available  in  two 
dimensions  but  it  should  be  noted  that  the  Stone -Weinst rass  theorem 
assures  the  ability  to  uniformly  approximate  a continuous  function  of 
2 variables  by  a bivariate  polynomial  with  arbitrary  accuracy  [2-16, 
2-17]  . 

Recently  Hou  [2-18]  and  Peyrovian  [2 -IP  ' have  employed 
spline  functions  to  image  restoration  problems  with  considerable 
success.  However,  in  both  of  these  works  the  knots  were  both 
fixed  and  uniformly  spaced  with  no  attempt  being  made  to  improvi  the 
approximation  by  adiusting  the  knot  placements.  In  one  dimension 
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the  existence  of  best  approximating  splines  with  both  fixed  and  free 
knots  have  been  shown  to  exist  [2-20"!  , and  that  the  approximation 
capabilities  of  splines  are  greatly  increased  with  the  allowance  of 
free  knots.  Recently  Schultz  has  developed  error  bounds  for  multi- 
variate spline  approximations,  in  both  and  norms  in  terms  of 
the  maximum  mesh  width  of  the  knots  and  the  moduli  of  continuity 
[2-21, 2-22]  . 

While  Schultz'  results  are  mainly  concerned  with  error  bounds 
and  not  the  best  approximating  properties  of  splines,  it  is  one  purpose 
of  this  dissertation  to  demonstrate  experimentally  that  a sensible 
placement  of  the  knots  adds  considerable  power  to  the  image  approxi- 
mating capability  of  bivariate  cubic  splines.  Furthermore,  since  the 
mapping  of  the  sampled  image  into  its  spline  coefficients  represents 
a nonenergy  conserving  transformation,  a bandwidth  reduction 
(manifested  as  a reduction  in  the  number  of  spline  coefficients  needed 
to  describe  the  image)  that  is  not  necessarily  at  the  expense  of  an 
increased  dynamic  range.  Thus  the  possibility  exists  of  more 
efficiently  quantizing  and  encoding  those  coefficients  than  those  of  the 
more  widely  used  orthogonal  transforms. 

2 . 3 Overview  of  the  Dissertation 

The  main  features  of  this  dissertation  are:  a)  the  adaptation  of 
a cont  i nuous -di  scret  e model  for  image  restoration  in  general  and  its 
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application  to  the  tomographic  imaging  system  in  particular;  b)  the 
use  of  the  eigenvalues  of  the  gram  matrix  in  predicting  the  tomo- 
graphic imaging  system  performance;  c)  the  adoption  of  an  approxi- 


mation theoretic  approach  to  relating  the  degrees  of  freedom  of  a 
sampled  image  to  its  original  unsampled  version  when  a system 
weighting  function  is  either  unknown  or  of  secondary  importance  to  the 
particular  image  processing  task;  and  d)  the  application  of  that 
approximation  theory  to  the  design  of  "smart  sensor"  imaging  arrays. 

The  following  is  a brief  summary  of  the  dissertation  on  a chapter 
by  chapter  basis. 

Chapter  1 serves  as  an  introduction  to  the  notion  of  image 
degrees  of  freedom  and  as  such  outlines  the  difficulties  involved. 
Models  for  imaging  systems  are  introduced  along  with  their  appropri- 
ate constraints  that  will  be  used  throughout  the  work. 

Chapter  2 serves  as  a review  of  the  research  efforts  towards 
characterizing  the  degrees  of  freedom  of  imaging  systems  with  a 
known  weighting  function.  The  recent  efforts  in  algorithm  develop- 
ment and  information  content  in  the  tomographic  imaging  system  are 
also  covered.  The  main  intent  is  to  review  known  systems  and  the 
approximation  theoretic  approach  will  be  reviewed  here  and  expanded 
upon  in  Chapter  5 after  it  has  been  properly  motivated. 

Chapter  3 contains  the  development  of  the  tools  necessary  for 
the  investigation  of  information  content  in  imaging  systems.  Here  the 


relation  between  the  gram  matrix  eigenvalues  and  system  degree  of 
freedom  is  developed  along  with  algorithms  for  estimating  f using  the 
gram  matrix  of  the  continuous -discrete  model.  This  is  done  for  both 
the  case  where  the  system  of  equations  is  invertible,  and  where  the 
system  is  singular  and  a constrained  least  squares  approach  is 
necessary.  The  separable  gramian  is  also  discussed  and  shown  to 
result  in  a significant  computational  reduction.  Finally,  the  relation 
of  the  gram  matrix  eigenvalues  to  the  continuous -continuous  model 
eigenvalues  is  explored  and  bounds  are  given  for  the  separable 
gramian. 

Chapter  4 is  the  application  of  the  results  of  Chapter  3 to  the 
tomographic  imaging  system.  A continuous -dis crete  model  along  with 
its  associated  gramian  is  developed.  Structure  in  the  gramian  is 
shown  to  exist  making  a linear  algebraic  solution  possible.  The 
degrees  of  freedom  of  the  tomographic  system  are  obtained  using 
the  gramian  and  are  shown  to  be  in  excellent  agreement  with  the 
general  type  of  system  tomographic  imaging  describes.  Numerical 
results  include  some  excellent  reconstructions  for  experimental 
projection  data. 

Chapter  5 is  devoted  to  the  information  content  of  a sampled 
image.  In  this  chapter  the  problem  is  shown  to  be  an  approximation 
problem  related  to  the  prior  methods  (where  the  approximating 


functions  are  taken  to  be  linear  combination  of  the  PSF  kernels  in  the 
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continuous -discrete  model).  However,  in  this  case  there  is  con- 
siderably more  flexibility.  The  concept  of  an  "epsilon  degrees  of 
freedom"  is  defined  and  related  to  sampled  images  using  splines 
with  free  knots. 

Chapter  6 is  devoted  to  the  experimental  results  of  the  "free"  or 
"variable  knot"  splines  of  Chapter  5. 

Chapter  7 contains  a summary  of  the  dissertation  along  with 
some  conclusions  and  possible  future  work. 

Appendix  A deals  with  the  relevant  computational  properties  of 
normalized  B-splines. 

Considering  the  dual  nature  of  this  work  a flow  diagram  of  the 
chapters  is  given  in  Fig.  2.  1.  It  should  be  useful  in  visualizing  the 
directions  taken  through  the  dissertation. 


k. 


Ch.  1 


Figure  2.1,  Organizat ion  of  the  Dissertation. 
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Chapter  3 

THE  CONTINUOUS-DISCRETE  MODEL  AND  THE  GRAMIAN 
3.  1 General  Development 

In  this  section  we  relate  the  DOF  of  the  imaging  system  (mani- 
fested as  dependent  data)  to  the  linear  dependence  amongst  the  PSF 
kernels  of  Eq.  (1-2).  In  doing  this  we  will  select  an  arbitrary  unit 
length  vector  a which  when  dotted  with  the  imaging  equation  yields 

<a,£>  = <a.  A*  h(?,  T))flMl)d^d7l>  + < a,n>  . (3-1) 

R 

If  the  first  term  on  the  right  is  much  smaller  than  the  second,  then  its 
contribution  can  be  considered  minimal.  Taking  this  with  the  fact  that 
'! Cl|1  = 1,  we  must  reduce  the  degree  of  independence  or  number  of 
DOF  by  one  for  every  independent  a such  that 

!<a  . |T  h(-.H)f(Ml)dPdTl  >|2  5 |!n||2  (3-2) 

R 

2 

where  ||  n|]  can  be  considered  the  sensor  noise  which  contributes  to 
the  definition  of  the  signal  to  noise  ratio  of  the  imaging  system.  We 
can  further  investigate  the  implications  of  equation  (3-2)  by  noting 
that  the  vector  inner  product  on  the  left  can  be  related  to  the 
continuous  inner  product  of  f with  h by 
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2 2 
|<a,  \\  h(§,H)f(5.‘n)dFdT>|  = |jj  f(?,Tl)<  h(f,Ti)(  a > d?dn  | 

R R 

and  applying  the  Schwartz  Inequality 

I JJ  fr,TU<h(?,Ua  >d^dn|2  sjy  |f(F,Ti)|2d-dTijj  |<h(Mi),a  >|2d-dT 

R R R 

= ||*(?.H)|£  a::[r]a 

Here  [pi  i9  the  gramian  of  the  PSF  kernels  and  is  defined  as 

[r]  = rf  h(?,71)h*(?,71)d!rdT 

■u 

R 

or 

Y.  = IT  h.(?,Tl)h:::(?,  TDdFdTi 

V)  J ' 1 J 

R 

which  is  Hermitian  having  real  eigenvalues  only.  Because  f(ctn)  can 
be  considered  bounded,  | f(F,  71)  | 5 E,  we  will  assume  in  conjunction 
with  equation  (3-2)  that 

EV[r]a  S ||  n |i2 

Expanding  [F]  into  its  eigenvector  decomposition  [f]  = [u1[a"|[Uj 
and  let  £.  = [u"|  a,  another  unit  vector,  we  have 

eVMI  * llall2  • 

Clearly  _3  *)  ^ is  minimized  by  the  allowing  of  t o be  the  eigen- 

vector  associated  with  the  smallest  eigenvalue  of  fT'1,  which  leads  to 


the  conclusion,  that  for  every  eigenvalue  X such  that: 


one  degree  of  freedom  is  lost. 

We  can  now  turn  our  attention  to  the  problem  of  obtaining  an 

A 

estimate  of  the  object  f,  and  demonstrate  the  role  played  by  the 
gramian.  In  doing  this  we  will  determine  a minimum  normed  least 
squares  estimator  with  the  constraint  that  the  norm  of  the  error  be 
equal  to  the  norm  of  the  noise  vector.  The  least  squares  approach  is 
taken  as  it  leads  to  a generalized  inverse  solution  in  the  discrete- 
discrete  model  which  Rao  £3-1,  3-2"]  has  shown  to  be  a minimum 
variance  unbiased  estimate  of  the  original  object,  f.  The  discrete- 
discrete  case  is  developed  first  and  serves  to  motivate  a minimum 
normed  least  squares  solution  to  the  continuous -discrete  problem. 
The  discrete -discrete  problem  is  nothing  more  than  a regression 
problem  wherein  the  image  vector  £ is  related  to  the  object  vector  f 
through  a blur  matrix  [h]  as  follows 

1 = [H]  f_  + n (3-3) 

Here  all  vectors  are  Mxl  column  vectors  and  [ill  is  an  MxN  matrix 
with  M^N  in  general,  n is  a noise  vector  whose  elements  will  be 
taken  to  be  identically  distributed  independent  random  variables  so 
t hat  H f ri  £ [3  - 2 1 . An  estimate  f can  be  obtained  from 


the  following: 


subject  to 


II  [Hit -4I|*  = ||  n |g 

which  results  in  the  minimization  of 

w(n  = iTf  + Y[([H]f-1);::([H]£-£)  + ||  n ll^} 

This  leads  to 

C [l]  + VrHT::[H]]£  = Y[h]::£ 

resulting  in 

i = [7  tl]+  [Hfa 

which  can  also  be  expressed  as  [3-3*1 

1 = [HJ  " [!]  + [H]  [H]  "J  *£  (3-4) 

The  Moore -Penrose  pseudo  inverse,  hf,  is  given  by  [3-3,  3-4"] 

[H]+  = lim  [H]*["-  [l]+  [H][Hf 

a 1 

and  taking  f_  = [ 11*1  ^ results  in  a minimum  normed  least  squares 

solution  to  Eq.  (3-3).  Also  note  that  [ H]+  involves  the  matrix  [h][h] 

which  is  a discrete  form  of  the  gramian. 

The  continuous-discrete  model  can  be  treated  in  a similar 

manner  by  stating  the  problem  as 

...  r>  2 

minimize  j f (A)dQ 


‘ -1 


(3-5) 


2 Z 

subject  to  11  Jh(6)f(0)d8 -g11^  = i'n!^ 


where  0=  (r,  Tl  ) for  notational  simplicity.  The  criterion,  UJ(f(6)),  to 
be  minimized  is  given  by 

ot(f( e )>  = J f2(0)de  + y f ||  Jhfeifteide-^  + |!  n ||2}  (3-6) 


Minimizing  tu( f ( 0 ) ) necessitates  the  use  of  the  calculus  of  variations 

A A A 

wherein  we  let  f(0)  = fQ(0)  + cA(e).  fQ( 6 > is  the  estimate  minimizing 
Eq.  (3-6)  and  A ( G ) is  any  function  of  6 . Thus  we  have  that: 


1)  u>(f(0))  = minimum  value 

e = 0 

A 

achieved  at  the  f ( 0 ) such  that 


2)  — uu  (f( 0))  = 0 

de  e=0 


ui(f(0))  * 0 


to  guarantee  a minimum.  Substituting  f(8)  = fg(0)  + G A ( 6 ) into  Eq. 
(3-6)  we  obtain 


UL)(f ( 0 ) ) = J ( f Q ( 0 ) +eA(e))2de+y^1£-2Y£TJh(0)(fo(6)  + eA(e#dQ 


+ yJJh'VeiMTKf  (0)  + eA(0)  ) • (yn)  + eA(TO)d6dTl  + y t|  n| \ 


(3-7) 
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I he  first  derivative  with  respect  to  e results  in 

If  ^(0))=  2J  (V0>  + sA'9))A(e)de  - 2y  £TJh(0)A(e)de 

+ 2yjh'  ni)(f0(Tl)  + eA(Tl))dT1  J h(G)A(0)d9 

^2 

7T  li,(*(e))  = 2Ja  ( 0 )d 9 + 2 y jj  h'(0)h(Tl)A(0)A{Tl)dedT1  :>  0 

os  ““ 


(3-8) 


If  y s o we  Will  have  a minimum-  Setting  Eq.  (3-8  ) equal  to  zero  and 
evaluating  at  e = 0 we  obtain  a minimum  at  fQ ( 0 ) satisfying 

JA(9)  {fo(0)  - Vh  : (9)^(0)  +Y_h  ' (0)J  h( T| )fQ ( Tj )d T1 } d 0 = 0 
which  will  be  satisfied  if 


?o(0)  + Yh  (e)Jh(Tl)f0(iri)dT1  = yh*(0)£ 

Notice  that  this  expression  is  of  the  algebraic  form 
(I  +YV)f0(Tfl  = Yir'(e)£ 

where 

^0(71)  " J h*(0)H(Tl)£o(Tl)dTl 

which  can  also  be  written 


(3-9) 


where 


VfQ(Tl)  = jK(e,  71)f0(T))d1 


K(0,  T1)  = h ::(0)h(Tl) 


1 

l 
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so  that  we  recognize  Eq.  (3-9)  as  a Fredholm  integral  equation  of  the 
second  kind  which  can  be  solved  by  iterated  kernels  giving 


OS 

fncn)  = £ (--yV)l6h  (e>£ 


oo  ^ 

X (-y^1  = (i  mi 

i = 0 


(3-10) 


if  (I+yV)  exist  s and  !!  yV  I!  £ 1 . 

To  obtain  an  explicit  form  for  the  solution  and  conditions  for 
convergence  note  that 


V f _ ( T) ) = I?  h=:  (e>h(  ?)h'  (P)h(T)fA (TDdTId- 

0 — — — U 


= f h*(6)[r]  H(T)  )f0(Tl)dTI 


implying  that  the  kernel  of  this  transformation  K (9,  7)  ) is  given  by 

k2  (an)  = h:(e)[r]h(Ti) 


and  that  V f(T])  can  be  written 


vnf0cn)  = J h:(9)[r  ]n‘1h(n)f()ni)di 

n „ 1 

and  that  K (9,V)  = h::;(0)[r  J " h(l)).  Thus 


fQ(0)  = Y \ ( -YV)  >h  (n»£ 

n-*oo  [ i =0  J 

= Yl'mjS  (-Y'7  h (0)rrl’'1h(n)h'(n)dT(j£ 

n-*t»  ' i =0 


ye)  = y h :(e>  rtr1  + v [r  l]'1  £ (3-in 

if  the  limit  exists.  However  since  rll  and  are  both  matrices  and 

frl  has  nonnegative  eigenvalues  only,  every  eigenvalue  of  ^il  + y r~" 
will  be  nonzero  for  y>  0 implying  that  [[  il  + y rrT  exis's.  Thus  we 
can  rearrange  the  right  side  of  eq.  (3-11)  and  obtain  the  estimate, 

f(?,  n)  = ye) 

n r1 

f(c,  T1)  = h (e,  D)  ~[l]+[r]  £ (3-12) 

Note  the  similarity  between  (3-11)  and  (3-9),  and  that  a Moore- 
Penrose  pseudo  inverse  will  be  obtained  by  letting  y-*»  so  that 

f(?,T)  = h>,T)[rf£ 

which  will  be 

f(*.  v)  = h " (-,  ,p)[r]‘1£ 

if  [r"]  is  nonsingular. 

A criterion  for  the  convergence  of  the  series  in  eq.  (3-10)  is 
t hat  IlyVll  = |t  y H V < 1 . This  has  physical  si gni  ficance  for  pas  s ive 
imaging  systems  since  the  output  image  energy  will  be  less  than  the 
input  image  energy  for  such  systems.  Thus 
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sllj  * II  f 


so  that 


llvf|L 


& 1 


resulting  in 


||vfM- 


MU  = sup  lifir 

2 MflUO 


$ 1 


At  this  point  it  is  interesting  to  investigate  the  relative  error  in  the 

A 

estimate  f(F,  71)  is  of  the  form 


f(F,71)  = <h:  (F,  71 ),  a> 


where 


a - [7  tn  + in] 


-l 


£ 


-l 


If  we  let  a'  = [u]  ^ so  that 


f(?.T])  = cp  ■ (?.  17),  a' 


where  the  cp.'  (?,  71)  are  orthogonal  and  c£"(F,  71)  = h_  ( F , 7))[u],  then  the 


relative  error  in  f(  = , 71)  will  be  given  by  • This  can  be 


!oMI 


accomplished  by  expanding  in  the  eigenspace  of  [r ] by  letting  [Ul  be 
such  that 

[u]*[r]  [Ul  = [Al 


where  [/\1  is  the  diagonal  matrix  of  eigenvalues  of  [r1  , and 
[Ul  ! [Ul  = [i]  . Thus 


2Q 


1 


££(  ' , T1 ) = [u]  ' h(e,  n ) 

and  the  solution  of  the  continuous-discrete  imaging  equation  will  be 
a solution  up  to  the  degrees  of  freedom  of  [P  "1  . Then  we  note  that 

T C£(-,  T%V,Tl)d*d11=  |T  [u]*h(-,  T))h  (',  P)[u]d?dn 

R R 

= [uf  [r]  [u] 

= [M. 

i.e.  the  set  fc3('r,T1)}  is  orthogonal.  Then  a'  is  as  follows 


a1  = tur^in+rnj1* 

= fy  [U]+  [u][A][uy::[ulj 

= [U]+  [U][A]]  £ 


If  [p]  is  nonsingular  the  pseudo  inverse  solution  is  obtained  by 


letting  y-*  ® so  that 


a'  = [a]"  [U] 


It  is  easily  verified  that  in  terms  of  C?(f , Tl) 


He,  n)  = <c2  ' (f,  n),  a'>  = h'l*,  7n[pf  £ (3-13) 

The  results  of  equation  (3-13)  indicate  an  ill-conditioning  problem 
in  which  an  upper  bound  on  the  relative  error  'I  5 a ' ' 1 / 1 a ' " can  bo 


found  t o be 


llSa'II 


In  general  then  we  can  obtain  our  estimates  for  f (?,  Tl)  in  the  eigenspace 
of  the  point  spread  functions  with  the  same  results,  and  can  generalize 
to  other  quadratic  constraints  by  letting  [c]  be  a nonsingular  constraint 
matrix  and  postulate  the  problem  as  follows. 


minimize  a *tc]a 

subject  to  Mj>-  [u]  [A  "I  a1^  = 0 . 


The  criterion  to  be  minimized  becomes 

a*[c]  a + y Ufi-tu]  [ A]  aj!2 

when  y is  the  usual  Lagrange  multiplier.  Performing  this  minimiza- 
tion results  in  a constrained  least  squares  solution  of 

a=  (Y'1[C]+  [aIV^AKuVs  (3-14a) 

and 

f(F,  n)  = a '[Ul  ' h(ff,  H)  . (3-14b) 

Proper  conditions  on  [C]  and  y result  in  the  pseudoinverse  reconstruc- 
tion and  are  developed  elsewhere  [3-5,  3 - 6 [1  . 

3.  2 Separable  Kernels 

In  image  processing  we  are  often  confronted  with  very  large 
matrices  that  make  implementation  of  the  algorithms  of  section  3.  1 
difficult  if  not  impossible.  If  it  is  noted  that  if  a sampled  image  data 
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point  is  obtained  at  each  coordinate  tuple  (x  , y. ) and  that  there  are  N 

11  x 

and  samples  in  the  x and  y directions  respectively,  then  the  image 

vector  g is  an  N *N  xl  column  vector  and  frl  will  be  a N N xN  N 
x y *-J  xyxy 

matrix.  Thus  for  images  of  even  a moderate  dimension,  say  1 00 

_ 4 4 

by  100,  IT]  will  be  a matrix  whose  dimension  is  10  x 10  containing 

g 

10  elements.  Clearly  some  structure  in  [Flmust  be  found  that 
allows  a simplification. 

The  structure  investigated  in  this  section  will  be  that  of  a 
separable  kernel.  By  a separable  kernel  it  is  meant  that  h . ( =■ . "P, ) can 
be  written  as  the  product  of  two  functions  as  follows: 

h.(F,Tl)  = h(1)(Z)h(2)(T})  Vi  = 1,2 N *N  (3-15) 

ill  x y 

This  subject  is  of  interest  becaus  e it  often  occurs  in  practice  and  that 

separability  allows  a significant  computation  reduction  over  the 

general  case. 

From  (3-15)  the  kernel  ve ct or  h(c,  T|)  of  (1-2)  becomes 

h(p,p,)  = hfl  V)®  h(2)(n) 

where  <g>  denotes  the  Kroneck.er,  or  direct  product  between  h*  * ' and 


h(p,m  = 


h(11)(?).h(2)(Tl) 


>(?)- h(2)(Tl) 


ii1)(S).h(Z)(1,) 

x N «N  xl 

L-  x y 

The  computational  savings  provided  by  kernel  separability  lies  in  the 
fact  that  [r]  can  be  shown  to  be  the  Kronecker  product  of  the  two 
gramians  [r]'  and  [Tj  ; as  follows 

[r]  = [r](1)®  [r](2) 

where 


[rf] 

= J h(1)(F)h(1)  (F)d§ 

(3-1 6a) 

[rf 

= J h(2)(Tl)h(2)  (Tt)d-n 

(3-1 6b) 

To  show  this  consider, 

[r]  = JT  h(5,  T|)h ' (f , Tl)dFdTl 
R 

which  for  separable  kernels  becomes 


[r]  = If  {h(1)(?)0  h(2)(T1)]{  h(1)(n®h(2)(H)fdPdTl 
R 

From  Bellman  [3-7l  we  know  that 

{h(1)(n®h(2%)}:  = h(1\n  ®h,2\n) 


thus 


[ri=  rr  {h(1  fe) ® h(2)(T!)}  {h(1)'(F)®  h(2)  (t )}  d^dr 

R 

But  by  the  algebra  of  Kroneckers  [3-7] 

[hfl)(c)®  h(Z)(71)]  [h(1)  (F)®h(2)  (Tl)}  = [h(1  V)h(1  ' (':i}©[h2(TJh(2)  (1] 
so  that 

[r]  = r h(1)(-)h(1)  (^d^®'  h,2'(Tl)h(21  (P)dP 

R_  r„ 

? Tl 

[rl  = [r](1)®  [r](2)  (3-17) 

Clearly  it  is  a significantly  simpler  task  to  diagnoalize  two  matrices 

of  dimension  N and  N each  than  one  N N xN  N matrix.  By  way  of 
x y x y x y 

illustration,  since  [l~"f  ' and  [r*^  ' are  Hermitiar^  they  are  diagonaliz- 
able  as  follows 

[aT*1’  = [u](1)  [r](1)[u](1) 

[a  12)  = [u](2)’[^](2)[uj2) 

then  from  the  algebra  of  Kroneckers  [3-7*1 

[a]=  [a](1 ’®[Af ' = [uf ' [r](1)[u](1)®[u](2)  [rl(2)[u](21 

= [uf  ’ ® [uf  ’ [r](  1 ' ® [r](2 ] [u]' 1 1 ® [ul(2 1 

and  finally  from  Bellman  [3-7]  we  have  that 


[[r](1)®  [rfY1  = [rf]  ® [r  f]~ 

Consequently,  the  coefficients  for  the  estimate  of  f become 

a = [A](1)  \vf)  ® [A f]  (3.18a 

f(Mi)  =/[r](1)  h{l\n®  [rf]  \{Z)cr\)  o-ist 

For  the  case  where  [r]  is  singular  and  the  constrained  least  squares 
solutions  are  used,  we  have  the  following 

i = [y-1[c]+[a](1)®  [a](2)H'1[a](1)[u](1)®  [a](2)[u](2)£ 

If  we  let  [c]  be  separable  and  nonsingular  as  is  often  the  case 
(take  for  example  [c"]  = [l]),  then 

a = CY"1[C](1)®[c](2)+[[Al(1)®[A](2)f]'1[A](1)[u](1)®[AT(2)[u](2)£ 


= [c](1)  ®[c];2)  [Y ‘1&]  + [c](2)  *[a](1)“®[c](2)  *[a](2)V 

• [a](1)[u](1)®[a](2)[u](2)£ 


•1  . 2 


1 , 2 


l = a;:;[u](1)  ®[U](2) 

Thus  we  see  that  the  constrained  least  squares  solution  to  a problem 
with  separable  f Y "]  and  [C^]  matrices  involves  the  diagonalizat  ion  and 
inversion  of  far  smaller  dimensioned  matrices.  The  discussion  of 
another  attractive  structure  is  deferred  until  the  section  dealing  with 
the  tomography  gramian. 


■ 


3.3  Gramian  Eigenvalues  and  System  Eigenvalues 

In  the  preceding  section  the  degrees  of  freedom  of  the  continuous- 
discrete  model  of  a linear  imaging  system  has  been  identified  with 
the  number  of  effectively  nonzero  eigenvalues  of  its  gram-matrix. 

This  prompts  the  question  of  exactly  what  is  the  relationship  between 
the  gramian  eigenvalues  and  the  original  continuous  system  eigen- 
values, for  it  would  certainly  seem  reasonable  that  they  exhibit  a 
similarity. 

In  this  section  we  will  apply  a result  of  Keller  ^3  — 8*1  that  will 
allow  us  to  relate  the  eigenvalues  of  [r^j  to  the  singular  values  of  a 
four  continuous  variable  kernel  h(x,  y;£,  n)  which  is  not  necessarily 
hermitian.  T o do  this  we  form  the  auxiliary  kernels  h'(x,  y;F,  7])  and 
h"(x,  y;1,  T|)  and  formulate  an  outer  product  expansion  for  h(x,  y;P,  T)). 
Under  the  original  assumption,  namely,  that  h(x,y;",7])  is  square 
integrable,  i.e., 

AT  |h(x,  y;*,  V)\  dxdyd-dT  < °° 

R 

and  defining  the  adjoint  kernel  h "(e,  lj;x,  y)  as 

h ' (ff,  T!;x,  y ) = h(x,  y;?,  71) 
then  the  auxiliary  kernels  are: 

h'(x,  y;F,  H)  = JT  h(x,y;n.e)h  (o.P;r,T)dnde 

= J*f  h(x,  y;p , 0)h( P,  Tt ;o . 0)dpd9  (3-19) 
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h"(x,  y;P,Tl)  = JJ  h:'(x,  y;p,  9)h(o,  6;?,  T))d0de 


and 


= JJ  h(p,  6;x,  y)h(p,  6;?,  T) )do <3 0 

(3-20) 

|0.cp.(x,  y)  = 

JJ  h'(x,  y;F,  T| )cp^. ( f . TDd^dTi 

(3-21) 

)i.i|r  (x,y)  = 

JJ  h'(x,  y^.TlW^.TDdfdn 

(3-22) 

1 

|u.  |2C(l(x,  y) 

= JJ  h(x,  y;F,T))t(F,Tl)dFdTl 

(3-23) 

|p.  |2^(x,  y) 

= yy  h:i(x,  y;£,  1))cp.(f,  TlldFdTl 

(3-24) 

Then  h(x,  y ; P , T] ) admits  of  the  following  expansion 


h(x,  y;*,  Tl)  = l2«C^(x,  yHy  ^ , n) 


(3-25) 


If  we  sample  h'(x,  y;§,  Tl)  at  the  coordinate  tuples  x^' ' = (x. , y. ) and 


x^1' ' = (x.,y.),  i - 1,...,N  N , j = 1, 
~ J J x y 


N N , it  should  be  clear 

x y 


t hat 


h'(xU\x'j1)  = y.. 

— — ii 


Let  be  a two  dimensional  quadrature  rule  approximating 

fff(x,y)dxdy  given  by: 

R N N N N 

X y x y 

ji  - y y f(x.y.)w y^  w . = i,  n = n *n 

..  4-f  *— { I j n 4—*  4^  ii  X y 

N * N 1 = 1 1 = 1 J J 1 1 

r <r  ** 


x y 


Jl*  f(x,y)dxdy 
R 
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Sampling  f(x,  y)  on  a cartesian  product  grid  we  can  separate  w into 
a product  of  two  weights  wjX',  wfy'  [3-9]  where  these  are  the  respec- 
tive quadrature  weights  for  one  dimensional  integrals  along  the  x 
and  y directions.  The  quadrafure  rule  then  becomes 

N 


Ji  - ^ f(x.,  y.)w'.X'w(.y' 

yn  2-«  1 J 1 


N xN 
x y i=l  j=l 


To  simplify  the  quadrature  notation  let  x = (x,  y),  0 = ( .-  , "p)  then  the 
samples  become: 


,(i) 

(x.,y.) 
t J 

i = 1 , 2,  . . . 

, N N 
x y 

xli]cR2 

(j) 

(5..TU 
J J 

j = 3,2,  ... 

, N N 
x y 

l(j)eR2 

,«>  = 

w!x).w!y) 

j j 

j = 1.2.  ... 

, N N 
x y 

Then  the  quadrature  approximation  to  (3-21)  becomes 

-(i)'  - £ 

e0» 


(3-26) 


or  in  vector  form 

uk®k  - 

where  [w]  is  a diagonal  matrix  of  quadrature  weights.  If  we  let  \ 

k 

t h 

be  the  k eigenvalue  of  [r  1 ^w"1  and  _U  its  associated  eigenvector, 
t hen 


considering  those  quadrature  rules  with  weights  greater  than  zero, 


and  such  that 


lim  J f = JT  f(x,  y)dxdy 

N N x N 

x x y 

N -»0B 

y 


Keller  has  shown  that  for  every  eigenvalue  X^  of  [r][w]  there  exists 
an  eigenvalue  4^  in  (3-21)  such  that 


where 


I 

lXk"Ukl  * Car+°(1)]2  max  | e ( x ) | (3-27) 

X 


|e(x)]  = | |J  h'(x;?,n)cpk(f,Tl)dFdll-  j h'm  | 

R N N 

x y 


and  is  the  area  of  R.  By  letting  R equal  the  unit  rectangle  and 
employing  a rectangular  integration  rule  with  evenly  spaced  samples 


such  that  w 


(j) 


1 


N .N 


, then  the  eigenvalues  of  [r  ] converge  to  a 


x y 

constant,  (dependent  on  the  number  of  samples),  times  some  eigen- 
value in  eq.  (3-21), 

A 

Since  the  error  in  the  estimate  f in  eq,  (3-6)  is  upper  bounded  by 

X 


the  condition  number  of  [r]  , |-maX  | , the  actual  value  of  the  cons- 
train 

tant  is  unimportant.  What  is  important  is  that  the  shape  of  the  eigen- 
value map  of  [r*1  rw]  will  be  in  good  agreement  with  the  spectrum  of 
the  kernel  in  the  continuous -continuous  model.  These  results  serve 
to  put  on  a firm  footing  what  should  be  our  intuitive  feeling  towards 


the  image  restoration  problem,  namely,  that  the  difficulty  in  restor- 
ation is  not  so  much  a function  of  the  sampling  method  used  (for  any 

one  consistent  with  a quadrature  rule  such  that  j f -*  I*  fdx  could  be 

n n J — 

used),  but  is  innate  to  the  original  continuous -cont inuous  formulation 
of  the  image  restoration  problem. 

3*  ^ Eigenvalue  Error  Bounds  for  Separable  Kernels 

For  separable  kernels  actual  numerical  bounds  on  the  error 

l\‘Ukl  can  obtained.  Wielandt  [3-10]  investigated  and  developed 

bounds  for  the  errors  in  estimating  the  eigenvalues  of  a hermitian 

operator  in  one  dimension  using  various  quadrature  rules.  Since  the 

two  dimensional  separable  problem  reduces  to  two  one  dimensional 

problems,  it  is  relatively  straightforward  to  extend  those  results  to 

the  two  dimensional  separable  problem.  In  the  previous  section  the 

eigenvalues  of  LI  1 [w]  were  investigated  where  w = — i — for  all  i 

i N N * 

a x y 

I his  presents  no  difficulties  in  obtaining  f in  Eq.  (3-5)  since  we 

simply  scale  both  sides  by  — ^ . The  results  of  the  previous 

x y 

section  are  also  valid  for  other  quadrature  rules  where  neither  the 
sampling  subintervals  nor  the  quadrature  weights  are  constant.  To 
deal  with  this  we  must  multiply  both  sides  of  Fq.  (3-13)  bv  the 
appropriate  quadrature  weight  matrix  rwl  and  obtain 

rw]£  = [w][r]fula  (3-28) 
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If  we  restrict  our  rules  to  those  with  the  property 


max  w. 

— : < M as  N N -*  « 

min  w.  x y 

i 

then  the  invert ibilit y of  (3-13)  is  dependent  on  and  not  jw"1.  Now 
[w]  [r]  is  similar  to  M[w]  which  is  itself  similar  to  the  hermitian 

i i 

matrix  [w]‘‘[r]  [w ~\i  implying  that  these  matrices  all  have  the  same 
eigenvalues.  Thus  we  can  analyze  the  separable  gramian  eigenvalues 
using  Wielandt's  results  applied  to  frl  ^w'  . 

As  before,  we  form  the  auxiliary  kernel  h'(x,  y;~ , T1 ) as  follows 


h'(x,  y;c',  Tl)  = h(x,  y;p,  c)h(r,  T;p , P)dpd« 

R 

= j*p  h^  \x,  p )h'  \y,  0)h  1 \ f , p)h*2  '(T,  0)dpd0 

R 

= |T  h*  \x,p)h(  '(T.o'dp"  h'^'(y,  6) h * ^ \ T) , G )d © 

R R 

X y 

Then 

(1)'  (?.)' 

h'( x,  y\~ , Tl)  h'  (x.  *)h'  (y,  T)  (3-29) 


For  a separable  kernel  Eq.  (3-21)  becomes 


M. 


I1  u^cpf1  \x)co'2)(y)  = J h'^'fx,  <r)cfl1)(,r)dtrJ'  hl2l'(y,  T)ci2'(T)d 

J J R J "r 

x y 


(3-30) 


whe  re 


(1)  (1)  f , (1)'  (1) 

u.  Cp:  (x)  = j h (x,  *)cp  (=r)d: 

R 1 


t t 


(3-31) 


1 
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(2)  (2)  » (2)'  (2) 
u:  cp.  (y)  = J h (y,  T )c2.  (Tl)dl 
R J 


I he  appropriate  quadrature  rules  with  x and  y sampled  on  a cartesian 
produce  grid  yield 

X-  U,  = IT]  rw>  U.  (3-33) 

1 L 1 

\^2\/2^  - t(2)  (2) 

Uj  = IT  J [wj  U.  (3  34) 

which  implies 

(D  (2>  _ (i)  (2)  V2'  (2) 

A : a.  - u.  qo  jj.  [ri  ® ri  . [wl  ® :w  .u  ®u 
1 J ^ J — i — j 

The  desired  result  is  to  bound 

(2)  (1)  (2.), 

X.  X.  -u.  u. 

t J ' J 

To  do  this  consider  that 

x,1\!2'-u,1,u!21  = n!1|-u!1|,o!2)-u!2|i  + ufl’u!2,-u|2'i 

1 J 1 J 1 > j i i i i 

(2)  (1 ) (1), 

+ u.  (X.  - u.  ) 

J i i 

(2)  ( 1 ) 

Since  u 2=  0 Vj  and  u.  ^ 0 Vi,  taking  the  triangle  inequality  implies 


Since  we  can  consider  'hat  J.‘]  si  and  u[2*  s 1 we  obtain 

t l 
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r 


] (i)  (2)  (i)  (2)  (i)  (i)  (2)  (2)  (2)  (2)  . (1)  (1) 

|X.  -u.  u.  I s | X.  -u.  ||X.  -u.  | + |X.  -u  |+|X.  -u.  I 

(3-35) 

(l)i 

Now  Wielandt's  result  applied  to  h (x,  ?)  for  rectangular  integration 

with  N points  'rom  To,  1 X equally  spaced  such  that  w = — , for  all 

x i N -1 

(1)'  X 
i is  such  that  if  h (x,  F)  is  a hermitian  Kernel  such  that 

1^  (x,  -)|  SL1  ' for  all  [O,  1 1 then  for  every  X*  * an  eigenvalue 

of  f there  exists  a u*  lan  eigenvalue  of  (3-27),such  tnat 

x 1 

, ,(1)  (1),  1.08L(1) 

1 \ "ui  I * “n 

X 

1 he  results  for  [T~f  1 are  identical.  Applying  this  to  (3-35)  we  have 
for  rectangular  integration  with  and  N evenly  spaced  points  in  the 
x and  y directions  that 


t , (1 ),  (2)  (1)  (2),  . (1.08)2  . (1).  (2)  1.08Lm  1.08L(2' 

X.  X.  -u.  u.  s — — L L 4 — — 4 — 

t J i J 1 N N N N 

x y x y 


(3-36) 


~~  h(1)'(x.  n s L(l1  Vr€  r0,  1 ] 
4x  * > ’ J 


~ h(?>  (y,  T| ) | £ I.'2'  \'Te  [0,  1 


(1 )'  (2)' 

If  the  kernels  h (x,  r)  ar.d  h (y,  T)  are  Iwice  differentiable  such 


h(1  ''(x,  ff)  si.'1'  VrcfO.  1 
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1 


3'  J2>« 

— j h <y.T) 

3y“ 


*l,2)  Vre  fo,  1 ] 


we  can  apply  Wielandt's  result  for  Simpson's  rule  to  (3-35)  and 
obtain 

i, (ll,(2|  111  (2),  . |.75|2l'‘’l'2’  .75l'‘’  .75L(2' 

lxi  Xj  -ui  ui  Is — ~i — ~t  * r 4 t 

J ‘ (N  -1)  (N  -1)  (N  -1)  (N  -1) 

X y X y 


(3-37) 


For  Gaussian  quadrature  over  a cartesian  product  grid  with  N 
and  N points  along  the  x and  y directions  respectively,  if  the 
auxiliary  kernels  h'  1 ' (x,  - ) and  h<2'  (y,  T>  ) possess  continuous  pth 


partials 


h,1,',X.O  and  A-  h|2>',y,n. 


,Pl 


P2 


?X*X  3y‘ 

then  for  N ^ p.  and  N 2 P-,  using  Wielandt's  results  for  Gaussian 
xi  y Z 

quadrature  it  can  be  shown 

(1 ), (2)  (1)  (2) 


|\:  \!  V u ’I  * W— -r) 

1 i j i l 1 \N  -1/ 


f dl  f 1 

/ d2  f2 

9P‘  .(I)'  e) 

VN  -1/ 

X 

In  - 1 ) max 

y x,  r 

i/.  ” ,x-  ’ 

• max 

y.  n 

d.  Pi 
.!  1 \r  1 

+ 4 1 — , ) max 


,N  -1/ 


d2  ^2 


4 4(n~T1 

y 


X,  ? 

max 

y.*n 


5y 

>Pi 


?P2  h'2>'l  T, 

— h ly.  Tl> 


3x 


Pi 


P2 


9yP2 


h(1)\x.’> 


h,2),(y.  Ti) 


whe 
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I 1.92  if  Pl,P2  *6 

drd2  { 

(4.0  if  pj , p2  i 7 

These  results  are  summarized  in  Table  3.  1 for  each  of  the 
aforementioned  quadrature  rules  in  both  one  and  two  dimensions. 

3.  5 Conclusions 

The  concept  of  the  gramian  eigenvalues  in  determining  the 
degrees  of  freedom  of  the  imaging  system  have  been  developed.  The 
gramian  eigenvalues  have  been  related  to  the  system  singular  values 
by  way  of  auxiliary  kernels  and  have  been  shown  to  be  closely 
related  to  those  singular  values.  Since  for  every  kernel  these 
singular  values  tend  to  zero, the  numerical  instability  in  digital 
image  restoration  arises  from  the  original  continuous -cont  inuous 
model  for  imaging  systems.  For  imaging  system  kernels  that  are 
separable,  actual  numerical  bounds  for  the  distance  between  the 
gramain  eigenvalues  and  the  square  of  the  system  singular  values 
were  found.  This  should  be  quite  useful  in  obtaining  upper  bounds 
for  the  number  of  effectively  independent  samples  one  can  obtain  for 
separable  imaging  systems. 
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Chapter  4 

PROJECTION  IMAGING  AND  THE  GR AMIAN  * 

4.  1 Continuous -Discrete  Model 

We  are  now  in  a position  to  apply  these  analytical  methods 
obtained  through  the  continuous -discrete  imaging  model  to  the  trans- 
axial  tomographic  projection  imaging  system.  In  the  ideal  case  the 
output  of  such  a system,  as  modeled  in  Figure  4,1,  is  a projection, 
p(r,  6),  related  to  the  original  object,  f(p,  71),  as  follows; 

p(r,0)  = JJ  f(f , T| )h( 0,  r ;£,  Tl)d cdTl . (4—1  i 

R 

In  any  real  system  finite  detector  width,  beam  spread,  scatter- 
ing and  other  phenomena  cause  the  line  projection  to  be  distorted. 
Thus  it  is  appropriate  to  consider  the  distorted  projection,  p(r,  0)  to 
be  related  to  f(f,  Tj  ) through  the  general  blur  h(9,  r;p,  Tl).  In  the 
important  case  where  the  blur  is  space  invariant  along  r,  and 
independent  of  0,  the  projection  becomes 

P(r,  6)  = JJ  f(?,  H)br  cos  P + Tjsin  P - rld^dT  . (4-2a) 

R 


This  and  portions  of  the  proceeding  chapter  are  summarized  in 
McCaughey,  D.  and  H.  C.  Andrews,  "Degrees  of  Freedom  in 
Projection  Imaging,  " I E E E Transactions  on  Acoustics,  Speech 
and  Signal  Processing,  vol.  25,  no.  1,  February,  1977. 
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where  b(?  cos  0+  Tlsin  0 - T1 ) is  b(r)  evaluated  along  the  line 


('  cos  9 + T)  sin  0 - r).  For  an  infinitely  narrow  blur  we  have  the 
ideal  imaging  system,  and  the  PSF  kernel  would  be  a line  mass  such 
that 

b(flr;?,Tl)  = 6(‘rcos0+Tlsin0-r) 

along  the  line  cos  0 + T|  sin  0 - r).  Hereafter  R will  be  taken  to  be 
the  unit  circle. 

For  this  ideal  imaging  system  the  line  mass  projection  becomes 

PL(r,0)  = ff  f(?,  Tl)6(f  cos  0+  I)  sin  0 - r)drdn  (4-2b) 

R 

and  it  can  be  shown  that  o(r,  0)  = p (r,  0)  'b(r).  If  the  Fourier  trans- 
form  of  o(r,  9)  is  denoted  by  Jfp(r,  0)},  then 

J[p(r,0)}  = J{pL  (r,  0)]jfb(r)}  (4-3) 

whe  re 

J[ p^r,  0)}  = J p(r , 0 )e  2n^Udr. 

- 00 

From  the  projection  slice  theorem  [4-1*1  we  can  relate  the  Fourier 
transform  of  the  projection  to  the  central  section  of  two-dimensional 
Fourier  transform  of  f(  - , Ti)  as  follows: 

j[p,  (r, pi]  = fj  re.P)..-*2™'’  ‘os  ,5'1  llsln 
R 

Denot  ing  3 {b(r  )j  as  R(u)  and  letting  g(u,c)  = ,7[p(r,ctj  we  have 


eiu.ei  = BiuijJ  c°s  e+  ”si"  SW.  ,4- 

R 

Finally  letting  g . be  g(u,  0)  evaluated  at  the  coordinate  tuple 

K,  l 

(u^,  A.),  we  have  the  following  continuous-discrete  characterization 
of  the  projection  imaging  process: 


gk,  i ■ B<V  II  e‘  + ”Sin  V d-dTI  (4-5) 


R 


i = 1,2 M 

0 

k = 1,2 

where  M and  are  the  number  of  samples  taken  around  0 and 

along  r,  respectively.  In  the  notation  of  the  previous  section,  then, 

the  total  number  of  samples  N must  equal  MM  = N. 

6 r 

By  lexicographically  ordering  the  indices  (i,k)  to  form  the 
image  vector  jj,  the  above  imaging  system  has  a gramian  consisting 
of  the  following  blocked  matrix 


tn  = 


tr]'1’1' 

[r](1>2) 

Iff-1' 

[r](2, 2) 

n(l.  Mo) 


tr]1 


,(M0,  1) 


mp) 


, . r_-i  (i , m) 

where  each  |_r  J is  a symmetric  M X M matrix  whose  k,  j l 


ent  ry  is  given  by 


- 


-i2TTuw("COS  ft+T  sin&  ) -fi^rm  ("cos0  +Tsin0  ) 
(i  ml  pp  J K 1 * l m m 

vjk.T,  = b<vb  <v  JT e 


R 


■ d "dT . 


Letting 

5 = r cos  <t> 

T]  = r sin  i> 
dxdy  = rdrd<5 

1 2tt  -i2TTu,  r cos(^-0.)  i2rru  r cos(<5-9  1 
(j  ml  ...  p p ^ k i l m , 

vli  ?/  = B(uk'B  '(ujl)J  j e e rdrd«S 


'(k,  J l 1 


0 0 


Letting  til  = <5-9.,  d^  = d 4,  -9.  s 2rr  - 9.  then 

1 2tt-9. 

•'•■'I  1 


U,  m)  m 

v(k.(1  =B,uk,B 


1 2tt-9  -j2mi,  r cos(f)  i2rru  r cos(ili+9.  - 9 ) 

k e'  1 1 m rdrd) 

e. 


Temporarily  let  A0  = 0.-Pm 


it-  ?:,  = B(«k)B:,(«()j  j 


1 2rr-9i  -j^Ttu  r cosft)  j2rru  r [cos(i|f  )cos (A0) 


'(k,  A) 


0 -0. 


-sin(f)ain(&9)]d|i 


d;>rdr 


= B(uR)B 


1 2n-9. 

!!'V  ‘ 

1 J0  -0. 


-j2nr  [cos  (ilf  )(u  cos  (A  9)-u,l  - si  nf'i'lu  s inf  A0)”! 

1 d^rdr 


Define: 


e 

c 


cos(A0)  - uk 


z 


e - u s i n A 9 
s A 


e = z cos(a),  e = z sinfal 
c s 

a = tan 


t 
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from  equation  (4-5)  provide  a structure  for  [ 1*1  which  facilitates  its 

generation,  diagonalization,  and  inversion.  If  we  note  that  0 and 

i 

0^  influence  [r]  only  through  the  cosine  of  their  difference,  by  samp- 
ling 0 uniformly  over  [o,  2tt)  and  u independent  of  0,  and  noting  that 

the  cosine  property  of  cos(0.-0  ) is  even  symmetric  about  M„,  rF"1 

l m 0r 

will  have  the  following  form  which  is  even  symmetric  circulant 
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: - 


only  in  the  indices  i,  m: 


tr]  = 


M(1’ 1 ' 

[r](1,2) 

[rf 1 - 3 ' 

[r] 

1 • 2 ) 

[r](1, 1] 

(rl11'2’... 

m(1>4) 

* j 

[r] 

[r](1’2) 

[rf’3) 

tr]'1,41... 

tr]'1-2' 

tr] 

(4-7) 


From  this  it  is  seen  that  each  row  of  blocks  is  a right  circular  shift 
of  the  preceding  row.  Thus  for  simplicity  we  can  let  [r]^’  m'  = 

[r  ] hereon.  In  the  special  case  where  each  [r]^  is  itself  a circu- 
lant,  [r]  becomes  a block  circulant  matrix.  Hunt  [4-2]  has  shown 
that  block  circulant  matrices  are  diagonalized  by  an  orthogonal 


matrix  that  is  the  kronecker  product  of  two  matrices  each  of  which 


represents  a discrete  Fourier  transform.  However,  for  our  pro- 
jection imaging  case  [t]^  is  only  symmetric  and  not  circulant  so  a 
more  general  result  than  Hunt's  is  required  if  we  are  to  effect  a 
computational  reduction  in  diagonalizing  [r]  . 

In  generalizing  Hunt's  result  to  include  the  pseudo-circulant2 
type  matrix  [1  ] we  are  confronted  with,  we  will  consider  a matrix 
[?],  whose  elements  are  taken  from  a complex  vector  space  V. 

I his  will  allow  a far  more  general  result  applicable  to  circulant 


2 

I he  term  pseudo-circulant  is  introduced  here  to  emphasize  that 
[rlis  neither  circulant  nor  block  circulant,  but  simply  has  two  of 
its  four  indices  operating  as  in  circulants. 
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matrices  of  rather  arbitrary  elements.  By  a complex  vector  space 
is  meant  a vector  space  whose  source  of  scalars  is  the  field  of  com- 
plex numbers,  C-.  Then  [f]  is  given  by 

Y(l)  ...  y(Me-l)_ 

Y(0)  ...  y(M  -2) 

Y<0) 

— 

where  Y( j ) are  elements  of  V and  may  be  scalars,  vectors,  matrices, 
or  tensors . 

Letting  UJ(k)  be  one  of  M roots  of  unity  the  closure  of  the  com- 

U 

plex  vector  space  V implies  that 


[r]  = 


Y(0) 


Y(Mp-l) 


Y(l) 


uu(k)y(i)€V  V Y(H6V,  tu(k)GJ 


Following  Hunt,  we  have  in  vector  form 


and 


uu(k)  = 


J2nk/Mg 


j2rrki /M( 


] = — — [u)(0)  ! 01(1)  ! ..  . | UJ(Mp-l  . 


(4-8) 


From  the  properties  of  the  complex  exponentials  the  matrix  [u^  "]  is 
orthogonal  and 


1 
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! 


K>  r-Mj]  M r%6] . 

The  matrix  [Ay]  is  a diagonal  matrix  whose  diagonal  elements  are 
elements  of  the  complex  vector  space  V which  may  in  fact  be  ma- 
trices. The  case  of  interest  here  is  when  the  vectors  y(i)  are  taken 
to  be  matrices  comprised  of  complex  numbers.  The  set  of  all  such 


matrices  clearly  forms  a complex  vector  space,  and  the  above  results 
apply.  Using  Kronecker  products  [4-3]  it  follows  that  the  matrix 

[r]  can  be  reduced  to  block  diagonal  form  by  the  matrix  [u^  1 $$ 

6 

[i]  and 


(4-9) 

where  [il  and  [a]^'  are  X matrices,  [i]  being  a diagonal, 

[A](°  being  full  and  being  defined  by  equation  (4-8).  The  block 

diagonal  matrix  on  the  left  side  of  (4-9)  contains  diagonal  matrices 
that  are  symmetric  if  [F](i'  is  symmetric  Vi  and  thus  diagonali zable . 
In  addition  the  orthogonal  matrix  reducing  the  left  side  of  (4-9)  to 
diagonal  form  is  itself  block  diagonal,  the  blocks  being  the  orthogonal 
matrices  diagonalizing  [a]('\  We  thus  have  the  result  that  the 
complete  diagonalizat ion  of  [r]  becomes 


1 


94 


£u]' 


i(M* 


[u] 


(2); 


M 


o 


o 


[ufMe): 


[A] 


(1) 


o 


tA]'2’ 

O ^'mq) 


|Cu] 


(i) 


rut 


(2) 


o 


(4-10) 


o 


,(Mp) 


where  fu"^  ' ^ [U  ^ for  all  i,  [d"^'  being  a diagonal 

matrix  of  complex  scalar  elements. 

If  Vi,  [r  is  a circulant  matrix.  Hunt's  result  for  block 
circulant  matrices  can  be  obtained  as  a special  case  of  (4-10)  by 
noting  that  the  set  of  all  circulant  matrices  whose  elements  are 
complex  numbers  comprises  a complex  vector  space. 

The  results  of  equation  (4-10)  can  be  summarized  as  requiring 
a general  discrete  Fourier  transform  to  reduce  the  block  diagonal 
form  followed  by  individual  subblock  diagonalizations  provided  by  a 
singular  value  decomposition  (SVD)  routine  developed  by  Golub  r4-4  . 
By  not  performing  the  multiplications  involving  zero,  a significant 
computational  reduction  is  achieved  in  diagonalizing  the  transaxial 
projection  imaging  gramian.  Furthermore,  the  nature  of  cos(  - ) 

over  [0,  2n)  implies  a further  computational  reduction  in  that  only 


5S 


an  even 


the  first  M^/2  + l blocks  of  [P  ] need  be  calculated  for  M 
number  since  the  first  row  becomes 

[r](1)[r](2)...[r](M0/2)[r](Me/2+1)[r](Me/2)  . r-^2) 

I hus  once  i = 1 » 2,  ....  M./2  +1 , is  determined,  [p]  is  known  in 

its  entirety.  This  also  applies  to  the  reduction  from  block  diagonal 
form  to  diagonal  form  since  the  discrete  Fourier  transform  is 
conjugate  symmetric  about  the  folding  frequency,  and  only  rA"l<l) 
i = 1 M0/2  +1  need  actually  be  diagonalized  by  the  SVD. 

lo  illustrate  the  significance  of  the*  above  simplifications  con- 
sider the  case  where  64  samples  are  taken  for  M and  512  samples 
for  Mp.  fp]  is  then  a 32,  768  dimensioned  matrix  containing  over  10 
elements.  Without  the  above  computat  ional  savings  it  would  be  un- 
feasible to  even  calculate  [r]  , let  alone  diagonalize  it.  However, 
by  calculating  only  the  first  M /2  4l  blocks  of  [p  ] , performing  the 
Fourier  transform  indicated  in  f4 . 9 > and  employing  the  SVD  algorithm 
of  Golub  [4-4"!,  this  particular  [r']  was  calculated,  diagonalized  and 
inverted  in  under  2 hour  epu  time  on  a PDP-10  computer. 

I he?  gramian  developed  in  the  above  transaxial  tomographic  form 
has  been  based  on  the  even  sampling  of  u (or  r)  and  resulting  in 
MyXMg  data  values.  While  this  even  sampling  of  1 and  has 
resulted  in  a significant  computational  reduction,  it  does  present  some 
difficulty.  In  discussing  this  and  for  tin  rcmiindi  r of  the  chapter  w. 


consider  the  projection  imaging  system  to  be  ideal;  namely,  p(r,  0)  = 


0^  (r,  R)  and  B(u)  = 1 in  equation  (4-4).  Looking  at  (4-5)  i'  can  be 

seen  that  if  u,  = 0,  for  some  k,  then  the  k-th  row  of  r"’1  ' will  be  the 
k 

same  for  each  i = 1,  . . . , Mf.  Therefore  the  rank  of  [r]  must  be  no 
greater  than  M f(  M ^ - 1 ) + 1 and  consequently  the  gramian  is  singular. 
Physically  the  u^  = 0 term  represent  the  d.  c.  or  average  value  of 
fp.Tl)  over  R,  which  of  course  is  rotation  invariant,  and  thus  no  new 
information  is  gained  concerning  the  d.  c.  value  after  the  first 
projection  is  obtained.  Naturally  this  will  result  if  a discrete 
Fourier  transform  is  applied  to  the  projection  data  to  obtain  g and 
the  u = 0 term  is  retained  in  each  projection.  Because  of  this  inher- 
ent singularity  in  [r  ' , the  estimate  of  f(lr,'T')  mus1  necessarily  result 
from  either  a pseudo-inverse  or  a constrained  least  squares  solution 
as  developed  in  equation  (3-14).  For  a minimum  norm  estimate  of 

A 

f ( F , T| ) we  would  set  the  constraint  matrix  such  that  r_C  1 = rl")  . 

Utilizing  a discrete  Fourier  transform  on  the  projection  data  to 
obtain  we  have 


E 


(1 ) 


(2) 


El  = 


(Mb) 

- 

where  g * is  an  M xl  column  vector  whose  k-th  element,  g,  '',  is  the 


k-th  harmonic  of  the  discrete  Fourier  transform  of  the  i-th  projecti  n. 
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I tilizing  this  result,  the  following  matrix  representation  of  the 
reconstruction  algorithm  is  obtained. 


= if  . hh ) ru ® fi]i 


tY'1[n+tA]2r1M 


[u]' 
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(4-11) 


rt  is  the  above  equation  which  is  utilized  as  the  reconstruction  algor- 
ithm for  the  images  in  subsequent  sections  of  this  chapter,  and  a few 
comments  are  in  order. 

Clearly  the  matrix  f]  *j  must  be  diagonalized  and  a constrained 
least  squares  inverse  precalculated  only  when  the  geometry  is 
changed.  Recalling  that 


fr,  T)  = h ( - , T)a 


and  that  the  [[  UJ,,  1 ® [i]]  can  be  implemented  using  fast  Fourier 
techniques,  the  bulk  of  the  computations  in  calculating  a lies  in  the 
multiplication  by  the  const  rained  least  squares  matrix.  Since  this 
matrix  is  block  diagonal,  by  not  performing  the  calculations  involv- 
ing zeros,  if  should  be  evident  that  this  step  involves  the  same 
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number  of  operations  as  the  convolutional  portion  of  a convolution 

algorithm  for  the  same  projection  data  format.  However,  the 

convolution  algorithm  involves  only  a single  sum  in  the  back 

projection  step  [4-1*1  while  this  algorithm  involves  a sum  over  both 

the  u,  and  0.  indices  where  k = 1 , . . . , M and  i = 1 Thus 

k i r 0 

the  convolution  algorithm  should  be  more  efficient  computationally. 

One  advantage  of  the  cont inuous -dis crete  formulation  described 

herein  is  that  it  provides  as  an  output,  a function  which  can  be 

evaluated  at  the  points  of  interest  without  an  interpolation  step  in  the 

algorithm.  This  eliminates  one  source  of  error  while  providing 

flexibility  in  the  output  format.  Reconstructing  a 64  x 64  output 

array,  approximately  8 min  of  CPU  time  was  required  for  M = 32 

r 

and  Mp  = 128.  I his  is  somewhat  slower  than  the  corresponding  4 min 
for  the  convolutional  algorithm,  and  for  larger  arrays  the  computa- 
tional advantage  should  favor  the  convolution  algorithm. 

4.2  Kxper  imentally  Determined  Degrees  of  Freedom 

To  further  emphasize  the  usefulness  of  the  gramian  in  estimat- 
ing the  degrees  of  freedom  of  an  imaging  system,  some  computational 
examples  are  developed  below.  T he  resulting  [F  "Iwas  determined 
and  diagonalized  by  the  procedure  described  in  the  preceding  section 

for  nine  different  combinations  of  M and  M . Three  values  of  M 

r 0 r 

equaling  16,  32,  and  64  were  selected  and  for  each  of  these,  three 
values  of  M,  were  selected  as  2M  , 4M  , and  8M-.  I hus  the 

r r 


1 


resulting  gramians  ranged  in  dimension  from  512  x 512  up  to 

.52,  768  x 32,  768,  Table  (4.  1)  lists  the  appropr  iat  e combi  nat  ions . 

The  gramian  associated  with  the  nine  combinations  of  M and 

r 

M r listed  in  Table  (4,  1),  were  formed  and  the  procedure  for  reducing 
the  full  [r  ] to  block  diagonal  form  using  Fourier  methods  ( [iuM  /)  ® 
fl  i was  employed.  The  block  diagonal  form  was  further  reduced 
utilizing  the  SVD  routines  of  Golub  ^4-41  , yielding  a completely 
diagonal  TV^M^  x matrix  with  eigenvalues  on  the  diagonal.  The 

eigenvalues  so  computed  are  plotted  in  figures  4.2,  4.  3,  4.4  for 

= ^6,  32,  and  64  respectively.  In  each  figure  results  are  shown 
for  Mq  equaling  4M^ , 8M^.  Since  nonsingular  permutation 

matrices  can  be  included  in  the  eigenvector  matrices  of  [fl  , without 
loss  of  generality  the  eigenvalues  are  plotted  in  decreasing  order. 

Note  that  due  to  the  large  dynamic  ranges  involved,  a logarithmic 
scale  is  used  and  that  XMAX  is  always  unity.  This  is  taken  to  be  the 
case  since  [r  1 can  always  be  normalized  by  dividing  by 
Figures  4.2,  4.  3 and  4.  4 display  a charade  rist  ic  behavior  in  the 

3 

A parenthetical  comment  is  in  order  here  for  those  readers 
familiar  with  tomographic  scanners.  In  this  ch.apter  we  have 
gathered  Mr  sample  points  along  each  radius  from  (0,  1)  and 
have  gathered  project  ions  at  angles  from  (0,  2rr  ).  In 
practice  projection  data  is  often  taken  radially  from  (-1.  1) 
and  at  angles  from  (0,tt)  (see  fig. 4. 1).  Clearly  the  same 
region  is  covered,  the  only  difference  being  a permutation  of 
indices . 
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spectrum  of  [r"j  . In  each  plot,  after  an  initial  decline,  X^  remains 

relatively  constant  with  increasing  n until  a critical  point,  called 

N , is  reached;  whereupon  an  abrupt  drop  of  4 orders  of 
crit 

magnitude  in  X occurs.  Subsequently  X^  decreases  more  slowly 
with  increasing  n.  However,  the  eigenvalues  were  calculated  on  a 
DEC  PDP-10  digital  computer  in  single  precision  arithmetic  and 

since  the  SVD  algorithm  calculates  eigenvalues  in  decreasing  ui-der, 

-6 

little  confidence  should  be  attributed  to  values  of  X less  than  10 

Consequently  we  can  attribute  the  tail  of  the  eigenvalues  spectrum  to 

be  due  to  computer  roundoff  and  computational  noise.  This  coupled 

with  the  abrupt  drop  in  \ at  N will  allow  us  to  roughly  estimate 

crit 

1 he  degrees  of  freedom  in  the  projection  imaging  system  to  be  . 

Note  that  this  abrupt  drop  allows  us  to  set  a threshold  of  approxi- 

_4 

mately  \ = 10  and  achieve  a consistent  mid  point  of  this  step  drop. 

-4 

Table  4.  1 includes  the  number  of  eigenvalues  above  10 

This  behavior  of  relatively  constant  X out  to  a point  followed  by 
an  abrupt  drop  has  been  shown  to  be  characteristic  of  ideal  circularly 
symmetric  systems  that  arc  bandlimited  and  observed  over  a circle  of 
finite  radius  [4-5]  . In  considering  the  projection  slice  theorem,  we 
note  that  the  Fourier  transform  of  the  output  projection  is  the  central 
section  of  the  Fourier  transform  of  the  original  cross  section. 
Therefore  even  ideal  projection  imaging  (i.e.  no  blur)  is  "band- 

limited"  by  inclusion  of  only  M harmonics,  and  thus  this  eigenvalue 

r f.l 


behavior  should  come  as  no  surprise.  More  will  be  said  about  this 
bandlimiting  process  later. 


In  figures  4,2,  4.3,  and  4.4  it  is  evident  that  increasing  M 

from  to  4M^  results  in  a significant  increase,  but  not  a doubling, 

In  ‘scr{t*  This  behavior  is  born  out  in  table  4.  1.  Also  in  table  4.  1 

2 2 

the  value  tt  (M^-l)  /4  is  included  with  other  parameters  of  interest. 

. 2 2 

Note  that  the  number  tt  (M^-l)  /4  is  in  relatively  good  agreement 
with  the  value  at  which  N does  not  significantly  increase  with 
increasing  Mg,  and  that  N is  always  less  than  the  maximum 
possible  rank  of  [r].  This  agreement  is  related  to  the  bandlimiting 
mentioned  before  and  can  be  explained  by  noting  that  the  continuous- 
continuous  ideal  projection  imaging  system  described  in  section  4.  1 
when  band-limited  with  a one  sided  bandwidth  B and  constrained  to 
a circular  input  region  of  radius  R,  is  a circularly  symmetric 
ideal  imaging  system  that  is  bandlimited  and  observed  over  a finite 
area.  Gori  and  Guattari  [4-6]  have  shown  that  the  degrees  of 
freedom  of  such  a cont inuous -continuous  system  are  in  reasonable 
agreement  with  the  Shannon  number,  which  for  a circularly  symmetric 
system  of  radius  R and  bandwidth  B has  been  shown  to  tie  *t  R r . 

In  our  case  if  w->  consider  the  pupil  to  have  Fourier  harmonics  no 
higher  t han  (M  - 11/2  and  take  R equal  to  unity  we  have  the  following 


DOF' 


estimate  of  the  degrees  of  freedom  or  N 


2 , 2 
TT  (M  -1  ) 
r 


DOF 


Note  that  this  estimate  of  the  maximum  number  of  effectively  nonzero 
eigenvalues  of  [r]  is  a function  of  M only  and  not  M„.  Intuitively 
this  makes  sense  since  increasing  samples  in  0 more  densely  fills 
the  circle  of  nonzero  frequency  components  of  radius  (M  -D/2  in 
frequency  space  and  does  not  increase  its  possible  area.  This  is 
exactly  the  situation  in  the  continuous  case  for  the  circularly 
symmetric  system  where  the  accessible  Fourier  area  is  a function 
of  the  radial  bandwidth  B.  More  rigorously,  the  results  of  Section  3 
have  shown  that  for  sufficiently  large  M^,  the  rectangular  integration 
rule  implicit  in  even  sampling  provides  a Gramian  whose  spectrum 
approaches  that  of  the  continuous -continuous  model  of  the  ideal 
projection  imaging  system  whose  aforementioned  degrees  of  freedom 
are  independent  of  9 . 

If  M is  specified,  considering  that  N s:  M • M we  should 
r crit  r 

certainly  not  expect  an  increase  in  reconstruction  performance  if 
M,is  taken  much  greater  than  N,^,„/M  . This  results  in  an  estimate 
for  the  number  of  projections  once  is  specified  as  follows: 

2 (M  -l)2 

n r 

M £ — 

0 4 M 

r 


which  for  sufficiently  large  reduces  to 


4.  3 Experimentally  Determined  Projections 

The  previous  section  was  devoted  to  investigating  the  degrees  of 
freedom  structure  provided  by  the  eigenvalue  spectrum  of  the  gram- 
ian  of  the  imaging  system.  In  this  section  we  turn  our  attention  to  the 

A 

actual  reconstruction  algorithm  to  estimate  the  object  f(?,Tj). 

Shown  in  figures  4.5,  4.6,  and  4.7  are  perspective  plots  of  a 
reconstructed  phantom  for  equaling  16,  32,  and  64.  As  before 
multiple  values  of  M are  taken  for  each  M . The  phantom  consists 
of  two  circular  regions  of  density  of  0.  1 and  0.  05  superimposed  on  a 
circular  background  of  density  0.  1.  These  results  we  obtained  with 

the  algorithm  of  equation  Klllwhere  y is  such  that  the  dynamic  range 

- 1 2 4 

of  t he  element  s of  the  diagonal  mat  rix  [y  As 

expected  from  the  eigenvalue  plot  for  M =16  an  increase  from 

M = 2M  to  M = 4M  produces  an  improvement  in  the  reconst  ruct  ion 

0 r e r 

while  an  increase  from  M = 4M  to  M = KM  has  a far  lesser 

r r 

effect  as  shown  in  Fig.  4.5.  However,  for  the  cases  of  M equaling 
32  and  64  the  respective  increases  from  M = 2M  l o M = 4M^ 
produce  little  noticeable  improvement.  Ibis  behavior  is  not  predicted 
by  eigenvalue  plots  for  M equaling  32,  and  64,  nor  should  it  be. 
t hese  plots  in  reality  predict  an  upper  limit  to  the  improve  ment  in 


A 


the  reconstruction  is  a function  of  an  ini  reased  number  of 


projections  given  the  number  of  sample  points  per  projection.  For 
the  cases  of  equaling  32  and  64  the  full  eigenvalue  spectrum  of 
the  matrix  [r^  is  probably  not  necessary  to  achieve  a reasonable 


reconstruction  of  this  relatively  simple  image.  This  fact  is  some- 
what emphasized  by  the  results  of  Figure  4.8  where  for  the  case 
of  equaling  64  a very  reasonable  reconstruction  is  obtained  for 
M equaling  64.  In  analyzing  this  it  should  be  considered  that  in 

U 

these  experiments  was  specified  a priori  at  3 different  levels  and 
for  this  simple  phantom  equaling  32  and  64  probably  indicates 
radial  oversampling.  This  would  seem  to  indicate  that  the  radial 
sampling  rate  should  be  set  to  reflect  the  highest  spatial  frequency 

believed  present  in  the  data.  Thus  if  M is  determined  in  this 

r 

manner  a more  substantial  improvement  in  the  reconstruction  should 
be  evident  in  increasing  M.  from  2M  to  4M  , but  in  no  case  will  a 
comparable  improvement  be  exhibited  by  increasing  M to  4M^  to 

8M  . 

r 

Figure  4.9  shows  the  reconstruction  of  a monkey's  head  which 
153 

was  irradiated  by  a radioisotope  source.  Hence  = 32  and 

projections  were  taken  at  1 degree  increments  from  0 to  360  degrees. 

Subsets  of  these  projections  were  taken  to  obtain  equaling  2M  , 

0 r 

4M  , and  8M  . In  this  case  increasing  M„from  2M  to  4M  produces 
r r 0 r r v 


an  improvement  in  the  reconstruction  while  an  increase  from  4M  to 

r 
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8M ^ produces  an  improvement  of  a comparably  lesser  degree. 
Obtaining  the  projections  by  the  method  described  probably  introduced 
some  angular  error  and  the  improvement  of  Fig.  4.9c  over  Fig.4.bb 
is  due  more  to  these  errors  being  averaged  than  any  significant  in- 
crease in  information  content. 

4.  4 Conclus  ions 

This  chapter  has  attempted  to  present  a generalized  degree  of 
freedom  analysis  for  imaging  systems  by  utilizing  the  gramian  of 
the  point  spread  function  kernel.  The  gramian's  eigenvalue  spectrum 
provides  an  indication  of  the  DOF.  and  various  sensor  signal  to  noise 
ratios  and  computational  noise  considerations  become  useful  para- 
meters in  the  eigenvalue  space.  The  genernl  gramian  approach 
to  imaging  was  then  directed  to  the  specific  imaging  geometry  of 
transaxial  tomographic  projections.  The  associated  gramian  was 
shown  to  have  considerable  structure  allowing  computational  savings 
in  calculation  of  the  eigenvalue  spectrum.  These  computational 
savings  utilized  a fast  Fotirier  routine  to  reduce  the  gramian  to  a 
block  diagonal  form.  This  block  diagonal  form  was  then  further 
reduced  by  performing  a series  of  lower  dimensional  singular  value 
decompositions  resulting  in  diagonalization  procedures  to  handle 
extremely  large  sized  gramians.  This  linear  algebraic  approach 
also  provided  a useful  reconstruction  algorithm  for  obtaining 
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estimates  of  the  original  object  from  the  projection  data. 

Experimental  verification  of  the  analytic  work  was  then  developed 
by  computational  procedures  on  nine  different  sized  gramians  for 
different  projection  imaging  sampling  geometries.  The  computational 
eigenvalue  spectra  so  obtained  agreed  quite  well  with  the  theoretical 
predictions  and  also  demonstrated  some  additional  structure  in  the 
imaging  system  inherent  in  circular  imaging.  The  degrees  of  free- 
dom were  quite  easily  obtained  from  the  eigenvalues.  A critical 
number  of  eigenvalues  was  consistently  observed  and  in  retrospect 
this  number  agreed  quite  well  with  the  inherent  bandlimit  of  such 
systems  when  analyzed  from  a continuous -continuous  imaging  model. 

Finally,  the  reconstruction  algorithm,  which  inherently  makes 
use  of  the  gramian  data,  was  exercised  for  obtaining  estimates  of 
the  original  objects  from  their  projection  data.  Again  experimental 
consistency  provided  added  confirmation  of  the  results  from  the 
gramian  analysis.  Essentially  the  visual  quality  of  the  pictorial 
reconstructions  agreed  quite  well  with  the  predicted  behavior  based 
upon  the  degree  of  freedom  analysis. 
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Figure  4.  8,  Phantom  Recon 


Chapter  5 

THE  DEGREES  OF  FREEDOM  OF  SAMPLED  IMAGE 
5. 1 Introduction 

In  the  preceding  chapter  the  number  of  effectively  nonzero 
eigenvalues  of  the  gram  matrix  has  been  used  to  characterize  the 
degrees  of  freedom  of  a sampled  image  when  viewed  as  the  output  of 
a linear  imaging  system.  However,  difficulties  arise  if  the  image  is 
considered  as  an  entity  to  itself  as  in  the  coding  problem.  While  it 
is  reasonable  in  the  discussion  of  sampled  images  to  assume  that 
we  are  dealing  with  bandlimited  scenes  sampled  at  least  at  the 
respective  Nyquist  rates  along  the  x and  y directions  in  the  output 
palne,  we  will  thus  assume  that  every  image  is  the  output  of  a linear 
imaging  system,  namely  the  ideal  low  pass  spatial  filter  whose 
Fourier  transform  H(u,  v)  is  as  follows: 

H(u,v)  = Rect(2i-)Rect(^-) 

x y 

where 

( 1 

Rect(x)  = \ 

* 0 otherwise 

B and  B are  the  one-sided  bandwidths  in  the  x and  y directions, 
x y 

The  system  weighting  function  is  then  given  by 
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h(x,  y ; p,  T] ) = 2B  Sinc(2B  (x-?))2B  sinc(2B  (v-Tl)) 

X x y y 1 


where 


sincfx)  = 


sin(nx) 


While  this  assumption  is  certainly  reasonable  it  leads  to  trivial 
results  as  will  now  be  shown. 

We  can  consider  the  bandlimited  image  f(x,  y)  under  consideration 

to  be  related  to  some  other  non-bandlimited  image  f'(5,  Tl)  through  the 

ideal  bandpass  filter  h(x,  y;£,  T])  as  follows 
00 

f(x>  y)  = [f  2BX  sinc(2Bx(x-P)£By  sinc(2B  (y-Tl))f'(?,  T)  )d£dT)  (5-1) 
which  results  in  the  following  continuous -dis crete  representation 

eo 

f(Xi’V  = If  2Bx  sinc(2Bx(xi-?))2B  sinc(2B  (y. -T|))f'(?,  T))dr dT) 

-oo  * y j 

Note  that  the  kernel  of  Eq.  (5-1)  separates  implying  that  [r"|  = [rf1'® 
[r](  ^ defining 


[rfM  = [y1  ] 

1 » J 


[r],2,=  [V21  ] 

J m,  n J 


(!)  . (2) 

Y.  . and  y are  as  follows 
i,j  m,  n 


(1)  2 r00 

Yi,j=4BxJ  sinc(2Bx(x  -?:))sinc(2B  (x.-f  ))d^  i.  j = 1 

-00  J 


Ym!n  = 4By  J sinc(2B  (ym-H)feinc(2B  (yn-T))dT)  m.  n = 1,  2 N 
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Working  with  (r "f  ^ if  we  let  0 = x,  -5,  it  can  be  shown  that  v' ^ ^ is 

1 i.J 

given  by  the  following  convolution: 

00 

y!  • = ^B2  f sinc(2B  (x. -x.-0))sinc(2B  G)d0 

X J X 1 J X 


Note  that 


and 


= 2B  sinc(2B  (x.-x.)) 
x x i y 


Y."! 
i.  J 


2B 


if  x.  = x. 
i J 

x.  -X. 

,f  nr1  ~ *’ 2>  • •• 


r 

m.  n 


2B 


if  = y„ 

m n 


y -y„ 

. , m n 

ff  " ‘t  2,  , , ( 


2B 


Thus  we  see  that  sampling  at  the  respective  Nyquist  rates  along  x and 
y produces  a gramian  that  is  an  N^.  square  matrix  that  is  diagonal 
as  follows 


tr]  = 4bxb  • xN  ® 


X X 


xN 

y y 


M = 4BA[I^  -N  XN  -N 
x y x y 


where  [l]^^  *s  an  NxN  identity  matrix. 


This  is  a little  less  than  satisfactory  for  another  reason, 
namely  that  to  obtain  this  full  rank  diagonal  gram  matrix  we  must 
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sample  at  the  respective  Nyquist  rates  everywhere  along  x and  y. 
Since  images  are  rarely  stationary,  the  frequency  content  of  the 
image  is  not  constant  over  the  entire  region  upon  which  the  image  is 
defined  which  would  intuitively  suggest  that  we  are  oversampling 
unless  the  frequency  content  of  the  image  is  spatially  invariant. 

Thus  to  define  the  degrees  of  freedom  for  an  image  of  this  type 
another  approach  is  necessary.  One  such  approach  is  to  apply  the 
singular  value  decomposition  (SVD)  algorithm  [5-1]  to  the  sampled 
image  matrix  whereupon  the  number  of  degrees  of  freedom  can  be 
equated  with  the  number  of  effectively  non-zero  singular  values  with 
additional  parameters  needed  to  describe  the  singular  vectors. 

This  method  deals  with  a sampled  version  of  the  image  only  and 
in  such  a way  that  the  degrees  of  freedom  are  affected  by  the  sampling 
method  used  and  can  result  in  misleading  results.  This  is  readily 
apparent  by  considering  an  image  f(x,  y)  that  can  be  written  as  the 
product  of  two  functions  fj  (x)  and  f2(x)  as  f(x,  y)  = f j (x)-  f2  (y ).  If  this 

image  is  sampled  on  a Cartesian  grid  (x.,  y ),  i = 1,2 N,  k = 1,2 

• • • * N,  then  the  image  matrix  can  be  written  as  the  outer  product 
of  the  two  vectors  [f{  (x, ). . . fj  (x^]  and  [^i  >•  • • VV"1  and  will  be 
a rank  one  matrix  for  all  such  separable  images  f.  This  represents, 
at  most,  2N  degrees  of  freedom.  Again  the  point  to  be  made  is  that 
the  degrees  of  freedom  should  be  a characteristic  of  the  original 


image  and  reflected  in  the  sampled  image  only  by  our  inability  to 
collect  an  uncountably  infinite  number  of  samples  for  application  on 
a computer. 

This  brings  up  some  conceptual  difficulties  since  the  number  of 
degrees  of  freedom  of  a function  defined  on  a continuum  is  countably 
infinite  at  best,  viz.,  the  space  of  all  square  integrable  functions  on 
[0,  I"]  wherein  any  function  can  be  expressed  in  a sense  by  a 
countably  infinite  expansion  of  any  orthonormal  basis  functions  that 
span  the  space.  Thus  we  need  to  determine  in  what  sense  a class  of 
functions  defined  on  the  unit  rectangle,  <£  = [x,  y ; -1  Sx,y  si}  are 
finitely  representable.  It  is  to  this  end  to  which  approximation 
theory  is  directed. 

5.  2 Degrees  of  Freedom  and  Shannon's  Sampling  Formula 

Before  proceeding  further  a few  words  considering  bandlimited 
functions  are  in  order  since  Shannon's  sampling  theorem  provides  a 
method  to  relate  the  original  function  to  its  sampled  version  and,  to 
some  extent,  the  ability  to  finitely  represent  the  image  when  it  is 
available  only  over  the  rectangle  fx,  y:  |x|  s X/2,  |y|  s Y/2}  . In 
this  case  the  image  f(x,  y)  is  given  in  terms  of  its  samples  taken  at 
the  respective  Nyquist  rates  along  x and  y as  follows: 
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f(x,  y)  = lim 

N -*• 

x 

N -»» 

y 


N 


x 

E 

i = -N 


x 


)sinc(2Bx(* 


(5-2) 


Here  B and  B are  the  respective  bandwidths  in  the  x and  v 
x y r 7 

directions.  Equation  (5-2)  involves  infinite  sums  and  in  any  practical 
situation  an  image  f(x,  y)  is  taken  to  be  non-zero  only  over  the 
rectangle,  so  that  = B^X  and  Ny  = B^y  results  in  4B^(B  Y 
samples  of  f(x,  y)  for  application  in  (5-2).  Then  finite  term  Shannon 
interpolation  becomes  an  approximation  problem  with  the  L ^ error 
given  by: 


E 

lil>Ny 


’ 2B  1 

y 


2 2 


y 

outside 

The  results  of  an  earlier  portion  of  this  section  concerning  the 
gram  matrix  of  the  bandlimiting  operator  are  interesting  in  this 
context,  in  that  they  indicate  that  sampling  at  the  respective  Nyquist 
rates  along  the  x and  y directions  will  produce  a diagonal  gramian  of 
constant  entries  whose  dimension  will  be  4B  XB  Y • Increasing  the 
number  of  samples  beyond  this  point  will  force  the  gram  matrix  to 
contain  non-zero  off  diagonal  terms  implying  an  increased 
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which  is  dependent  upon  the  rate  at  which 


i 


A 


correlation  between  samples.  Since  this  gram  matrix  is  that  of  a 
bandlimited  operator  observed  over  a finite  area,  increasing  the 
number  of  samples  sufficiently  will  produce  a gramian  eigenvalue  map 
similar  to  that  of  the  cont inuous -cont inuous  operat or- -namely  one 
that  is  constant  out  t a point,  not  significantly  greater  than  4B^(B  Y. 
with  an  abrupt  drop  beyond  that  point.  Thus  we  can  consider  that 
there  could  be  approximately  4B^XB^Y  independent  samples  of  f(x,  y) 
in  R . 

In  reality  Shannon's  sampling  theorem  provides  a method  of 
bridging  the  gap  between  the  continuous  and  discrete  for  bandlimited 
images  by  providing  a reconstruction  method  involving  uniformly 
spaced  samples  over  R . However,  by  considering  the  problem  in  a 
more  general  approach  it  may  be  possible  to  define  other  adaptive 
and  possibly  more  efficient  approximation  methods. 

In  this  light  consider  that  the  spatial  frequency  content  of  most 
images  is  not  constant  over  the  entire  image  and  the  possibility  exists 
that  if  the  rectangle  /p  is  broken  up  into  2 sets  ^ and  ^ such  that 
R-  U R^  then  the  portions  of  f(x,y)  defined  over  £ and  f can  be 
considered  sections  of  a bandlimited  image  with  x and  y bandwidths 
\ B^  ' and  B^  \ B^  in  each  of  these  regions.  By  way  of 
example  consider  ^ = fx,  y:  -X/2  SxSO,  |y|<  Y/2],  /p  - 
(x,  y:  0 £ x £ X/2,  |y  | £ Y/2}  . Then  f(x,  y)  need  be  sampled  only  at 
the  respective  Nyquist  rates  in  each  of  these  regions. 
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For  and  P,  respectively  we  obtain  = 2B^  ^XB  Y and 

12  x y x y 

(2)  (2) 

= 2B^  XB^Y  so  that  the  approximate  number  of  independent 

samples  is  2XYB  (B*  ^+B^  ^).  Now  either  B^  ^ £ B or  B^  ^ sB  So 
y x x x x x x 

that  it  follows  that  the  number  of  independent  samples  is  less  than  or 
equal  to  4B^XB^Y.  This  heuristic  argument  is  presented  to  motivate 
the  question  of  why  should  an  image  be  sampled  over  its  entire 
domain  of  definition  at  a Nyquist  rate  dependent  on  a bandwidth  that 
may  be  relevant  to  only  a small  region.  A similar  idea  will  be 
presented  in  the  next  chapter  for  bicubic  splines. 

In  summary  then  the  eigenvalue  map  of  the  bandlimited  gram 
matrix  gives  an  upper  bound  to  the  number  of  independent  samples 
available  and  a subsectioning  of  the  image  represents  a low  order 
attempt  at  reducing  this  upper  bound --hope fully  at  an  acceptable  error. 
5.  3 The  Degrees  of  Freedom  Viewed  as  an  Approximation  Problem 
In  characterizing  the  degrees  of  freedom  of  an  image  as  an 
approximation  problem  we  are  confronted  with  two  questions,  namely: 

1)  to  what  extent  is  f(x,  y)  finitely  representable,  and  2)  the  determin- 
ation of  the  finite  representation  from  a sampled  version  of  f(x,  y). 

In  dealing  with  these  questions  we  will  take  f(x,  y)  to  be  an  ele  - 

ment  of  a metric  space  W with  metric  d^, . For  example,  this  could 
2 

be  the  space  L (8 ) where  (x,y:  -lsx,  y si]  , if  we  do  not  dis- 
tinguish functions  that  differ  only  on  a set  of  measure  zero.  Consider 
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L 


then  the  following  generalized  bivariate  polynomial  approximation 


scheme  J : 
N 


N N 
x 


f(x,  y)  = Y]  V]  a°J.(x)l(x)  =/  N = N • N 
M 1 J N * Y 

with  the  property  that  if  we  have  any  other  approximation  scheme  J 


where 


that 


N N 

x 


•k  = § aij2i(xUj(y) 


Vf’V  £ 


dW(f’o/N)  * 


Then  is  a best  approximating  generalized  polynomial  in  the  metric 

V We  also  require  that  ■+f(x,y)  as  N ->  co  so  that  for  every  e>0, 

there  exists  an  N(e,  i,  d^)  (in  general  dependent  on  the  set  { l]  and  the 

metric  d ) such  that  N N >N(e  , A,dw)  implies  that  d,ir(f,^°  ) <e. 
vv  x y w VV  W 

Thus  we  could  define  the  degrees  of  freedom  of  f at  level 
epsilon,  or  more  succinctly  the  epsilon  degrees  of  freedom,  DOF 
(dw>  e ),  as 

DOF(d  , e ) = inf  {N(e,£,d  ):  d (f,y°)<e}. 

W W 


Here  the  infimum  is  taken  over  all  approximating  functions  not  equal 
to  f to  avoid  the  trivial  case  of  DOF(d^,  e ) equaling  1,  (the  approxi- 
mating function  being  f itself). 


94 


In  general,  this  is  a difficult  problem  since  for  each  set  of 

functions  { ^ fx)  ^ (y).  . . ^ (x)  jr  (y)}  we  must  find  a best  approxi- 

1 X 1 y 

mating  bivariate  generalized  polynomial  (if  it  exists)  and  then 

determine  the  set  of  functions  minimizing  N N . Thus  we  are 

x y 

concerned  with  the  existence  of  a best  approximation  in  addition  to 
its  determination. 

However,  if  we  consider  a one  dimensional  case  where  N is 

taken  to  be  the  space  of  bandlimited  functions  observed  over  a 

finite  interval  of  the  real  line,  the  distance  d„r  being  the  L-,  metric 

Vv  ° 2 ’ 

the  result  is  known.  Here  for  every  c>0  the  inf  N(e,  l,d  ) was 

W 

found  by  Landau  and  Poliak  [5-2]  to  be  achieved  by  the  functions  ($5} 
which  are  related  to  the  prolate  spheroidal  waveforms,  as  first 
described  by  Slepian  and  Poliak  [5-3"].  While  this  is  known,  the 
determination  and  utilization  of  these  waveforms  for  large  space 
bandwidth  products  has  met  with  little  success.  Most  likely  then  to 
arrive  at  any  significant  results  in  two  dimensions,  we  will  have  to 
restrict  the  search  to  those  sets  of  functions  that  are  computationally 
feasible  and  possess  the  approximating  properties  required. 

In  the  context  of  a digital  computer  the  actual  functions  available 
are  those  that  can  be  obtained  by  finite  sums  and  products.  Thus 
ultimately  all  functions  must  be  reduced  to  those  that  can  either  be 
generated  by  recursions  or  by  truncated  and  shifted  polynomials.  So 
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A 


a restriction  of  the  search  of  functions  to  polynomial  splines  would 


seem  to  be  reasonable. 

Polynomial  splines  are  chosen  due  to  their  approximation 
properties  and  the  fact  that  they  possess  a basis  namely,  the 
normalized  B-spline  basis,  that  provides  a local  basis  property  thus 
allowing  a rapid  generation  with  the  matrices  involved  in  generating 
a B-spline  fit  to  a function  f being  well  conditioned.  With  DeBoor's 
algorithm  for  computations  using  normalized  B-splines  [5-4]  no 
difficulties  are  encountered  in  handling  multiple  order  knots.  Here- 


after we  will  consider  a spline  sk  N N of  the  order  k with  (the  degree 

’ x y 

being  equal  to  k-1)  N and  N knots  in  the  x and  y directions  respec- 
tively, to  be  of  the  following  form: 

N N 

x 


k,  N N 
x y 


^ SiJNi.k(I!X,Nj.k^ 


y)  = ff x,  y) 


where  N , (•  ) are  the  normalized  B splines  of  order  k satisfying  the 
i,  k 

following  recursion  relationship  over  the  knot  sets  { F^ , F^, . . . , ^ ^ 

x 

and  {t1|,  . * * * ' ^ [■* 1 • 

y 

X-F 


N.  , (F;  x) 
i,  k — 


e.  , "x 

i+k 


i+k-1  i 


F Ni,k-1(-‘,X)  + F.  -F.,,  Ni  + 1,  k-l(-'’X) 


i+k  ‘ i+1 


N. 


,(!;*)  = 
i, 1 


1 xe[VFi+i 


0 otherwise 


Ob 


and 


N 

2 ■ 1 Y“t?l FN  1 

1 = 1 X 

These  functions  are  discussed  in  Appendix  A. 

For  the  duration  of  this  chapter  we  will  consider  f(x,  y)e 

with  the  distance  being  the  L2  metric.  Thus  for  any  x and  y knots 

sets  { 5 j . . ■ f ^ } and  f Tl^  , . , . , Tl ^ 1 the  matrix  [s.  .1  of  coefficient  s 
x y -1 

minimizing  the  L.,  error  for  a particular  set  of  x and  y knots  is  given 

by  the  solution  to  the  following  matrix  representation  of  the  normal 

equations : 


Ts  1 Trl(y) 

JN  xN  L ij  JN  xN  L‘  xN 
x y x y x y 


= IT  ii  ;x)N^(Tl;y)f(x.  y)dxdy 


(5-3) 


where  N^(^_;x)  and  N^(T[;y)  are  respectively  the  Nx  and  N column 
vectors  of  the  normalized  B -spline  basic  functions  and 

[r],x)  = 


[r],y)  = tv!!"! 


,(y)> 


where 


Yij 


(X)  = f N.  , (i;x)N. 


! j,  k — 


;x)dx 


Y. 


rj  = f N ■ , (T[;y)N.  ,(T;y)dy 
’J  ^ ' » k j , k 

Thus  finding  the  best  approximation  to  f(x,  y)  over  the  x and  y knots 

sets  can  be  stated  as  follows: 


>7 


Minimize 


■[ 


|f>  |f(x,  y)  - Sk  N (x,  y)|2dxdy 
S ’ x y 


■]* 


over  all  possible  x and  y knot  vectors  ( F ^ § } , . . . , T]  } 

x Ny 

subject  to 


I s 1 Vi  = 1 , 


N 


hjl sl  vj  = 1. 


N 


This  is  nothing  more  than  nonlinear  minimization  over  the  possible 
knot  sets  in  <?,  the  solution  of  which  has  been  shown  to  exist  by  Rice 

[5-5]. 

Thus  specifying  an  error  e , we  can  find  the  epsilon  degrees  of 
freedom  by  a sequence  of  minimizations  decreasing  the  number  of 
knots  in  the  x and  y directions  until  we  reach  a point  at  which  the 
error  will  be  exceeded  with  any  further  restriction  in  the  number  of 
knots. 

For  this  to  make  sense  we  must  be  assured  that  for  every  e>0 

there  exists  N and  N such  that 
x y 


f(X,y)  ' Sk  N N (x’y,^2  < e 
’ x y 


and  that  a minimum  over  the  knot  sets  exists.  Addressing  this  c 
vergence  problem  Schultz  [5-6]  has  shown  for  k = 4 that 


2 


lf,x,y)  ' SN  N (X’y)^  Sin)  il  4 f(X’y)||  + 2 2 f(x’y* 

x y " 1 9x  "2  dx  3y 

liHLl 


where  o = max(max(5,  max(T).  .-11.)}.  Thus  taking  N £—  , 

i i+l  i . j+l  J x p 

2 J 

N a — we  have  as  N , N -*  ® such  that  p -♦  0 the  bicubic  spline 
y - x y 

0 

approximation  will  converge  to  f(x,  y)  in  an  sense  over  the  unit 
rectangle  $ . 

5.  4 Computational  Considerations  and  Conclusions 

Heretofore  the  concern  has  been  the  existence  of  a solution  to 
the  problem  of  finding  DOFfd^e  ) of  an  image  on  the  unit  rectangle. 
While  a solution  exists  in  the  form  of  a sequence  of  best  approxi- 
mating splines  with  a successively  decreasing  number  of  knots,  this 
requires  the  solution  of  a nonlinear  minimization  at  each  step.  For 
approximations  of  the  form  of  (5-3)  involving  even  a modest  number 
of  knots  this  is  computationally  infeasible.  Furthermore  since  we 
are  dealing  with  a computer,  the  right  side  of  (5-3)  must  be  cal- 
culated using  numerical  quadrature  and  this  involves  only  an 
approximation.  Thus  we  must  settle  for  a procedure  that  provides 
an  approximation  with  a good,  if  not  optimal,  knot  placement.  In 
the  next  chapter  some  numerical  results  involving  bicubic  spline 
approximations  to  simulated  and  actual  images  will  be  given.  It  will 
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be  shown  that  there  exist  relatively  simple  methods  for  placing  the 
knots  that  provide  much  better  approximations  than  a bicubic  spline 
with  uniformly  placed  knots. 


Chapter  6 


EXPERIMENTAL  RESULTS  FOR  SPLINE  APPROXIMATION 


6.  1 Introduction 

The  purpose  of  this  chapter  is  to  present  some  numerical 
results  concerning  the  "epsilon  degrees  of  freedom"  concepts  devel- 
oped in  the  preceding  chapter.  This  was  developed  as  an  approxi- 
mation problem  whose  solution  was  seen  to  involve  the  determination 
of  a sequence  of  best  approximating  (in  an  L^  sense)  B-splines  with 
variable  knots.  While  the  determination  at  each  step  of  such  a best 
approximating  spline  is  simply  a nonlinear  minimization  problem 
over  the  knots  defining  the  spline,  it  is  computationally  infeasible. 
Thus  we  must  follow  DeBoor  [6-l"l  and  settle  for  spline  approxima- 
tions with  good  if  not  optimal  knot  placements.  In  what  follows,  two 
easily  implemented  methods  for  placing  the  knots  will  be  given  that 
can  result  in  a significant  error  reduction  over  the  uniform  knot  case 
for  the  proper  class  of  images. 

The  results  in  this  chapter  will  be  developed  using  cubic  splines 


givij.g  the  following  fourth  order  spline  approximation  f(x,y) 

N N 

r 

(x.y)  (6-1) 


= it!  siiNi  4^x>Ni 

i=i  j=i  lJ  l'*  h 


A s 


4,  N , N 
x y 
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r 


where  ? and  T1  are  the  knot  vectors  in  the  x and  y directions  respec- 
tively. The  spline  coefficients,  S were  obtained  by  solving  the 
following  system  of  equations: 


N N 
x y 


f(xrym' 


g g 

for  l - 1, 2,  ....  N and  m = 1,  2,  . . . , N.  In  matrix  notation  this 


becomes 


tf(xrym^NxN  " ^Ni.4(^SCAxNy^Sij^NxxNy^Nj,4(^,ym'1'NyxN 

(6-2) 


where  [ 1 indicates  matrix  transpose.  To  simplify  notation  let 


[f(xt,ym)]  = [F] 


[NT'  = 

[N]®  = [N.i4®yml] 


equation  (6-2)  becomes 

[F]  = [n1-)T[S..][N]® 

Since  N >N  and  N >N  in  general  equation  (6-2)  cannot  be  solved 

x y 

exactly.  However  the  spline  coefficients  minimizing  the  normalized 
least  squares  error  e given  by 
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1=1  m=l 


£ £ 

£=1  m=l 


(6-3) 


can  be  obtained  by  taking  [S_]  to  be 


[S..]  = [[t#1  [N^W  [F][N]®  [[N-fi*  [N^’r1 


(6-4) 


6.  2 Knot  Placement  from  Projections 

In  this  section  we  investigated  the  possibility  of  placing  the 

knots  for  the  x and  y knot  vectors  from  the  projections  of  f(x,y)  along 

the  y and  x directions  respectively.  To  do  this  we  borrowed  an  idea 

s t 

suggested  by  DeBoor  [6-1],  wherein  he  suggested  placing  the  i+1  — 
knot  with  respect  to  the  i— - knot  for  a k—  order  spline  according  to 
the  following 


J‘i+1  |f,k,(*)|1/kd*  = 


constant 


where  f^(x)  indicates  the  k^  derivative  of  f(x).  To  determine  the  x 
knot  set  for  a bicubic  spline  fit  the  knots  were  such  that 


?. 


V 


+1 


4 1 

— 4 [ f(x,y)dy 
7)x  ' -1 


dx  = constant 


(f>-5) 
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I he  y knot  set  is  obtained  in  a similar  manner.  Note  that  if  f(x,y) 


can  be  written  as  f(x,  y)  = fj(x)«  t y ) then  equation  (6-4)  can  be 
written  as 


r i+l 


ctx 


*4  fl(x) 


•if  1 


fz(y)dy 


dx 


and  we  have  that 


const  ant 


i 


4 

4 


i 

4 


(x) 


dx 


constant 

7 « 

f2(y)dy 

j 1 c- 


It  is  reasonable  to  expect  that  this  method  of  placing  the  knots  would 

be  quite  effective  for  images  that  are  either  separable  or  exhibit  a 

high  degree  of  separability,  and  not  so  effective  for  images  that  do 

not . This  point  is  illustrated  for  the  following  experiments  involving 

two  images--one  an  analytic  image  consisting  of  a Gaussian  pulse 

2 

with  standard  deviation  n equaling  0.  1 and  the  other  consisting  of  an 
actual  image  of  an  armored  personnel  carrier  (APC).  For  each 
experiment  N was  taken  to  be  128  so  that  the  effects  of  the  quadrature 
error  implicit  in  the  residual  of  equation  (6-3)  and  the  computations  in 
equation  (6-4)  should  not  be  a significant  factor  in  evaluating  the 
results. 


The  results  for  a bicubic  spline  fit  to  the  Gaussian  pulse  with 


10  knots  placed  by  the  projection  method  are  shown  in  figure  6.  1. 
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Also  shown  are  the  results  obtained  with  10  and  20  knots  placed 
uniformly,  the  original  image,  and  the  knot  placements  for  10  knots 
placed  uniformly  and  by  the  projection  method.  The  corresponding 
mean  square  errors  are  tabulated  in  table  6.1.  Note  that  visually 
the  results  of  10  knots  placed  according  to  this  projection  method  are 
as  good  as  of  20  knots  placed  uniformly  and  are  far  better  than  those 
obtained  by  placing  10  knots  uniformly.  From  the  results  of  table  6.  1 
it  is  clear  that  placing  10  knots  in  the  x and  y directions  by  the 
projection  method  is  better  than  placing  them  uniformly  even  up  to 
the  situation  employing  20  knots. 

Figure  6.2  shows  the  results  of  a bicubic  spline  approximation 
to  an  actual  image  of  an  APC  for  40  uniform  knots  in  each  direction 
and  for  40  knots  placed  by  the  projection  method.  Also  shown  are 
the  corresponding  knot  placements  in  each  direction  for  both  cases. 

Figure  6.3  is  the  original.  Here  the  results  indicate  that  the 
projection  method  is  not  so  good  as  the  uniform  method.  This  is 
evident  in  table  6,2  where  the  projection  method  results  in  an  error 
that  is  an  order  of  magnitude  greater  than  that  obtained  in  the  uniform 
knot  case. 

While  these  results  show  that  the  proper  placement  of  the  knots 
can  result  in  a significant  error  improvement  in  the  approximation  of 
separable  images,  another  method  must  be  employed  for  those  images 
that  are  not  separable. 
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6.  3 Spline  Approximations  by  Subsectionin 


In  this  section  the  possibility  of  subsectioning  the  image  and 
using  different  knot  densities  in  each  of  the  subsections  is  investi- 
gated. This  method  might  provide  fruitful  results  when  one  considers 
an  error  bound  given  by  Schultz  ^6 - 2 for  cubic  splines.  Recalling 
that  the  error  is  given  by  the  norm,  difference  between 

the  function  and  its  approximation,  this  bound  is  given  by 
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From  the  discussion  of  Chapter  5 concerning  local  bandwidth  and 
Shannon  sampling  it  follows  that  if  the  image  derivative  energy  is 
large  only  over  a small  region,  then  using  a uniform  knot  bicubic 
spline  with  knot  mesh  width  equaling  p given  by  equation  (6-7)  should 
result  in  an  overly  good  approximation  of  the  image  in  those  regions 
where  the  image  derivative  energy  is  low.  Thus  we  should  tie  able 
to  obtain  reasonable  results  by  employing  a different  bicubic  spline 
with  uniformly  spaced  knots  in  each  subsection,  t lie  knot  density  in 


I 06 


L 


each  subsection  being  proportional  to  the  value  of 


f(x,  y) 

+ 

*2  *2  „ , 
2 2 y) 

+ 

■Ar  f'x.  y) 

2 

Sx  By 

2 

sy 

dx 

in  that  subsection. 

This  approach  was  taken  for  the  APC  of  figure  6.  3.  The  image 
was  partitioned  into  64  square  sections  and  the  derivative  energy  in 
each  section  was  calculated  by  central  differencing.  Figure  6.4 


shows  the  results  for  this  method  along  with  the  associated  knot 
density.  It  can  be  seen  that  this  associated  knot  density  is  highest  in 
the  region  containing  the  insignia  and  lowest  in  the  sandy  regions 
surrounding  the  vehicle.  The  bright  lines  at  the  subsection  boundaries 
indicate  a fourth  order  knot  at  the  boundary,  and  note  that  the  regions 
containing  the  minimum  number  of  knots  have  x and  y knot  vectors 
consisting  of  two  fourth  order  knots  each  at  the  subsection  boundary. 
As  noted  in  Appendix  A this  corresponds  to  a bivariate  cubic  poly- 
nomial approximation  in  each  such  subsection. 

Table  6.  3 lists  the  number  of  parameters  necessary  to  re- 
construct the  image,  the  reconstruction  error,  the  data  reduction 
ratio  for  the  uniform  knot  case,  the  protection  placement  case,  and 
the  subpartitioning  case.  Note  that  the  error  is  lowest  for  the  sub- 
partitioning case  and  that  it  provides  the  best  reconstruction. 

To  further  explore  the  subpartit  ioning  method  a series  of 
bicubic  approximations  involving  image  subsections  of  different  sizes 
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was  run  on  the  APC  image,  an  aerial  reconnaisance  image  and  an 
image  of  Los  Angeles  International  Airport  (LAX).  For  this  series 
the  image  dimension  N was  taken  to  be  256  and  three  subsection  sizes 
of  32  x 32,  16  x 16,  and  8x8  pixels  were  used  for  both  images.  For 
each  subpartition  size  three  knot  density  ranges  were  employed.  In 
all  cases  the  maximum  knot  density  is  taken  to  be  such  that  the 
matrices  of  normalized  B-splines  in  equation  (6-4)  are  nonsingular 
thus  resulting  in  the  image  being  interpolated  in  at  least  one  sub- 
section. The  lowest  knot  density  in  each  subpartition  sequence  was 
taken  to  be  that  corresponding  to  a fourth  order  knot  placed  at  each 
of  the  subregion  boundaries  for  the  x and  y knot  vectors  respectively. 
The  number  of  knots  was  then  increased  by  raising  the  minimum 
number  of  knots  employed.  7'he  results  for  the  APC  image  along 
with  the  associated  knot  densities  for  subpartitions  of  size  32  x 32, 
16x16  and  8x8  are  shown  in  figures  6.  5,  6.  6,  and  6.  7 respectively. 
Figure  6.8  is  the  original.  Here  the  fourth  order  knots  at  the  sub- 
partition boundaries  are  not  displayed  for  aesthetic  purposes.  Note 
that  all  of  the  approximations  are  quite  good  and  that  the  knot 
densities  are  quite  adaptive  for  each  subpartition  size.  1 he 
corresponding  error,  data  reduction  ratio,  and  number  of  para- 
meters necessary  for  the  bicubic  spline  approximation  are  listed  in 
table  6.4.  Note  that  the  error  for  the  32  x 32  case  corresponding  to 
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the  data  reduction  ratio  of  5.62:1  is  lower  than  the  error  for  the 
16  x 16  case  for  the  5.  3 1 : 1 data  reduction  ratio.  This  would  seem  to 
indicate  that  for  low  error  levels  the  32  x 32  partition  size  is  better 
than  the  16  x 16  sized  partition  which  itself  is  better  than  the  8x8 
case. 

This  same  sequence  was  obtained  for  the  reconnaisance  image 
and  the  results  along  with  the  corresponding  knot  density  patterns 
are  shown  in  figures  6.  9,  6.10,  and  6. 11  for  the  32  x32,  16  x16,  and 
8x8  cases  respectively.  Figure  6.  12  is  the  original.  Table  6.  5 
lists  the  relevant  errors  and  corresponding  parameters  as  for  the 
APC.  Note  again  that  the  best  reconstruction  for  the  32  x 32  case  is 
at  a lower  error  than  the  16  x 16  case,  and  at  a higher  corresponding 
data  reduction  ratio.  This  would  again  seem  to  indicate  that  for  low 
error  levels  the  32  x 32  sized  partition  is  better  than  the  16  x 16 
partition. 

The  sequence  for  the  image  of  LAX  is  contained  in  figures  6.  13, 
6.  14,  and  6.  15.  Figure  6.  16  is  the  original.  Note  here  that  while 
the  errors  are  relatively  high  the  results  are  quite  good  visually 
and  that  the  knots  serve  fairly  well  in  locating  the  areas  containing 
the  aircraft.  Note  also  that  the  best  results  occur  for  the  8x8 
subpartitioning  case.  This  is  most  likely  due  to  the  relatively  small 
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size  of  the  items  of  interest  - in  this  case  the  aircraft. 


6.  4 Summary  and  Conclusions 

In  this  chapter  the  attempt  has  been  to  demonstrate  the  utility 
of  variable  knot  splines  in  achieving  a data  reduction  and  quantifying 
the  degrees  of  freedom  of  sample  images  by  the  number  of  variable 
knot  bicubic  splines  necessary  to  approximate  a particular  image  at 
an  error  level  epsilon.  Considerable  success  was  achieved  for  an 
analytical  image  consisting  of  a Gaussian  pulse  with  o = 0.  1.  Here 
10  knots,  each  in  the  x and  y directions  were  placed  by  an  algorithm 
dependent  on  the  4 partials  of  the  projections  and  provided  an  error 
reduction  of  three  orders  of  magnitude  over  the  situation  where  each 
of  the  knot  sets  were  uniformly  spaced  in  the  x and  y directions.  The 
error  was  such  that  10  knots  placed  in  the  above  manner  provided  an 
error  lower  than  that  achieved  by  placing  20  knots  uniformly  in  the  x 
and  y directions.  This  method  was  not  very  successful  for  the  APC 
image  due  to  the  lack  of  any  separability  characteristics.  However 
by  subsectioning  the  image  and  employing  a different  spline  approxi- 
mation in  each  subsection  whose  knot  density  was  dependent  on  the 
derivative  energy  in  that  region  good  results  were  obtained.  A high 
degree  of  adaptability  was  in  evidence  through  the  knot  density 
patterns  with  acceptable  errors  being  obtained  at  reasonable  data 
reduction  ratios. 
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NUMBER  OF 
KNOTS 

PLACEMENT  MODE 

MEAN  SQUARE 
ERROR 

10 

uniform 

2.  158xl0'3 

10 

4th  differences 

1.29  X 10'6 

20 

uniform 

3.24  x 10'6 

Table  6.1.  Mean  Square  Error  for  Knot  Placement  on  Gaussian 
Pulse  a =0.1. 


NUMBER  OF 
KNOTS 

PLACEMENT  MODE 

MEAN  SQUARE 
ERROR 

40 

uniform 

6.413  x 10~3 

40 

4th  differences 

-2 

6.  308  x 10 

Table  6.2,  Mean  Square  Error  for  Knot  Placement  on  APC  Image. 


PLACEMENT 

MODE 

NUMBER  OF 
PARAMETERS 

DATA  REDUCTION 
RATIO 

MEAN 

SQUARE 

ERROR 

uniform 

1936 

8.46:1 

6.  4 1 3x 1 o" 3 

4th  differences 
on  projections 

2016 

8.  12;1 

6.  308xl0~2 

sub  part  it  ioni  ng 

1885 

8. 69:1 

5.  57Hx 1 o" 3 

1 1 1 


T able  6.  3 


Data  Reduction  and  Errors  for  a 128  x 128  Bicubic 
Spline  Reconstruction  of  the  APC  Image. 


Parameters  for 
Reconnaissance  Image 


Bi- Cubic  Spline  Fit  with  10  Knots 

in  X and  Y Direction  Determined  from  4th 

Partial  s of  X and  Y Projections. 


Bi- Cubic  Spline  Fit  with  10  Uniform 
Knots  in  X and  Y Directions. 


.Figure  6.  1 Results  For  A Bi- Cubic  Spline 

Fit  To  A Gaussian  Pulse  With  =.  1 
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Bi- Cubic  Spline  Fit  with  20  Uniform 
Knots  in  X and  Y Directions 


Original 


Figure  6.  1 (continued).  Results  For  A Bi-Cubic 

Spline  Fit  To  A Gaussian  Pulse  With  0^=1 


] 


40  Knots  Uniformly  Placed  in  X 
and  Y Directon 


40  Knots  in  X and  Y Directions 
placed  by  4th  Differences  on  X and  Y 
Projections 

Figure  6.2  Bi-Cubic  Spline  Approximation  For  APC 
With  40  Knots 
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Parameter  Reducti 


Parameter  Reduction 


Parameter  Reduction 


MSP  .25°; 


Bicubic  Spline  Reconstructions  and  A 
Densities  for  an  APC  Photograph  Usi 
Size  32  by  32. 


Parameter  Reduction 
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Parameter  Reduction 


ieure  6.  6 Bicubic  Spline  Reconstructions  and  Associated  Knot 

Densities  for  an  APC  Photograph  l sing  Subregions  of 
Size  16  bv  16. 


Parameter  Reduction  5.31:1 
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Parameter  Reduction  3.62:1 


MSE  = . 17% 


Figure  6.7  Bicubic  Spline  Reconstructions  and  Associated  Knot 

Densities  for  an  APC  Photograph  Using  Subregions  of 
Size  8 by  8. 
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Parameter  Reduction  2.56:1 


and  Associated  Knot 
Photograph  Using 


12H 


Parameter 


Parameter  Reduction  = 2.35:1 


Parameter  Reduction  = 1 


Bicubic  Spline  Reconstructions 
Densities  for  a Reconnaissance 
Subregions  of  Size  8 bv  8. 


and  Associated  Knot 
Photograph  Using 


HHBV  i 

= 3. 31: 

1 

MSE  = .77% 

Figure  6.12 


Original  256  by  256  Reconnaissance 
Image. 


no 


Parameter  Reduct  ion  3 . 4 1 : 1 


MSE  = 2.8% 


Parameter  Reduction  5.48:1 


MSE  = 1.9% 


Figure  6.1  3.  Bicubic  Spline  Reconstructions  and  Associated  Knot 
Densities  for  LAX  Photograph  Using  Subregions  of 
Size  32  x 32. 
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Parameter  Reduction  9.  lQ:l 
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Parameter  Reduction  5.  3 5 : 1 
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Parameter  Reduction  2.82:1 


Figure  6.14.  Bicubic  Spline  Reconstructions  and  Associated  Knot 
Densities  for  LAX  Photograph  Using  Subregions  of 
Size  1 6 x 1 6.  132 


Figure  Bicubic  Spline  Reconstructions  and  Associated  Knot 

Densities  for  LAX  Photograph  Using  Subregions  of 
Size  8x8. 
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Chapter  7 


SUMMARY,  CONCLUSIONS  AND  FUTURE  WORK 
7.  1 Summary  and  Conclusions 

The  dissertation  has  presented  a degrees  of  freedom  analysis  of 
images  and  imaging  systems  in  digital  image  processing.  The  analy- 
sis has  been  applied  to  image  restoration  when  the  image  is  viewed  as 
the  output  of  a linear  imaging  system,  the  problem  being  to  restore 
the  original  image  up  to  the  degrees  of  freedom  of  the  system,  once 
this  has  been  quantified;  and  to  sampled  images  directly  where  the 
desire  is  to  relate  the  sampled  image  to  the  original  unsampled  and 
unblurred  version.  In  reality  both  situations  are  approximation 
problems  since  the  restoration  of  an  image  through  the  continuous- 
discrete  model  involves  the  determination  of  an  estimate  for  the 
original  image  in  terms  of  a linear  combination  of  the  point  spread 
function  (PSF)  kernels  of  the  model,  while  in  the  latter  case 
the  researcher  is  free  to  find  an  approximation  over  a wide 
class  of  approximating  functions.  In  modeling  the  restoration 
problem  the  cont inuous -dis crete  model  was  adopted  since  it  most 
closely  represents  the  imaging- -data  gathering-reconstruction 
sequence,  eliminating  any  quadrature  approximation  to  the  super- 
position integral  implicit  in  the  linear  system  assumption.  Further- 
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more,  the  eigenvalues  of  the  gram  matrix  of  the  PSF  kernels  were 
shown  to  determine  to  what  extent  restoration  of  the  image  is 
possible. 

The  assumption  that  imaging  systems  are  linear  and  represent- 
able by  a two  dimensional  Fredholm  integral  equation  of  the  first 
kind  presents  difficulties  in  the  restoration  process  since  restoration 
is  no  more  than  an  attempt  to  numerically  invert  the  associated 
integral  equation.  In  this  dissertation  we  have  bounded  the  error 
between  the  eigenvalues  of  the  system  gram  matrix  and  the  singular 
values  of  the  original  continuous -continuous  model  kernel.  It  was 
shown  for  a wide  class  of  quadrature  rules  with  their  corresponding 
sampling  schemes  that  the  eigenvalues  of  the  gramian  converge  to 
the  square  of  the  singular  values  of  the  continuous -continuous 
kernel.  Since  this  singular  value  sequence  possesses  zero  as  its  only 
limit  point,  it  is  obvious  as  to  why  the  image  restoration  process 
becomes  ill-conditioned  for  large  data  ensembles. 

More  important  is  the  fact  that  if  the  eigenvalue  behavior  of 
the  original  kernel  is  known,  it  is  then  possible  to  predict  the  number 
of  effectively  independent  samples  that  can  be  obtained  from  an 
imaging  system. 

The  continuous  discrete  model  was  developed  for  imaging 
systems  with  separable  kernels  and  the  associated  gram  matrix  was 


shown  to  be  the  direct  product  of  two  smaller  gramians  with  a 


resultant  computational  reduction.  Actual  numerical  bounds  were 
obtained  for  the  difference  between  the  gramian  eigenvalues  and  the 
continuous -continuous  model  singular-values  for  several  quadrature 
rules. 

These  concepts  were  applied  to  the  tomographic  imaging  system 
with  excellent  results.  The  continuous -discrete  model  was  shown  to 
possess  a structure  that  made  diagonalization  of  the  tomographic 
gramian  possible,  and  a linear  algebraic  solution  feasible.  The 
tomographic  imaging  system  was  shown  to  be  a circularly  symmetric 
imaging  system  that  is  bandlimited  and  observed  over  a finite  radius 
in  the  output  plane  and  whose  degrees  of  freedom  and  eigenvalue 
behavior  are  known.  The  gramian  eigenvalue  behavior  was  exactly 
as  expected- -exhibiting  a rather  constant  behavior  out  to  a predicted 
point  followed  with  a large  drop.  The  reconstruction  algorithm 
developed  from  the  continuous -discrete  model  provided  excellent 
reconstructions  from  real  projection  data. 

In  dealing  with  the  image  itself  the  degrees  of  freedom  was 
approached  as  an  approximation  problem  where  the  degrees  of 
freedom  at  a level  epsilon  was  taken  to  be  the  minimum  number  of 
functions  needed  to  approximate  the  image  with  an  error  epsilon. 

Since  this  minimum  is  difficult  to  find,  the  functions  used  win  cubic 
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splines  with  variable  knots.  By  dividing  the  image  into  subregions  a 
significant  data  reduction  was  achieved  with  reasonable  errors.  It 
was  found  that  the  number  of  knots  and  thus  the  degrees  of  freedom 
was  higher  in  regions  of  higher  image  derivative  energy  than  in 
those  regions  where  the  image  was  relatively  constant.  It  was  also 
found  in  one  case  that  while  the  mean  squared  errors  in  this  scheme 
were  of  the  same  order  as  obtained  with  a uniform  knot  spline 
approximation  with  a comparable  number  of  parameters  this  method 
produced  reconstructions  with  far  more  visible  detail. 

Finally  it  should  be  said  that  in  effect  this  represents  an 
attempt  to  bridge  the  gap  between  the  continuous  domain  upon  which 
images  are  defined  and  the  discrete  grids  upon  which  they  are 
sampled  and  defined  for  analysis  by  digital  techniques. 

7.  2 Future  Work 

There  are  several  topics  concerning  the  spline  approximations 
presented  in  this  dissertation  that  deserve  further  study.  this  is  not 
to  say  that  there  are  no  other  interesting  questions  concerning  the 
gram  matrix,  for  its  application  to  Fredholm  integral  equations  of 
the  second  kind  is  certainly  one  example.  Nevertheless  the  topics 
discussed  here  will  mainly  be  concerned  with  spline  approximations 
as  an  image  coding  device. 

In  this  light  we  have  seen  that  a sampled  image  can  be  approx- 
imated with  a tolerable  error  by  a bicubic  spline  with  fewer 
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coefficients  than  pixels  in  the  original  image,  and  as  such  it  repre- 
sents a data  reduction  scheme.  However  this  is  a reduction  in  the 
number  of  computer  words  and  not  bits  necessary  to  represent  the 
image,  and  in  order  for  these  results  to  be  applicable  to  digital 
communication  systems  they  must  be  in  terms  of  bit  reductions.  lhis 
would  indicate  that  an  efficient  quantization  method  for  the  spline 
coefficients  must  be  developed.  The  conjecture  that  the  dynamic 
range  of  these  coefficients  might  not  be  significantly  greater  than  that 
of  the  original  image  makes  this  a particularly  interesting  topic  for 
further  study.  Furthermore  the  normalized  B-spline  basis  functions 
exhibit  a local  basis  property  which  might  transfer  some  of  the 
image  spatial  correlation  properties  to  the  spline  coefficients  thus 
making  a further  de -cor relat  ion  possible  with  some  post -processing. 

In  this  work  the  image  was  subsectioned  into  square  regions 
only.  In  each  of  these  regions  the  knot  density  was  quite  adaptive  to 
the  image  derivative  characteristics  and  it  would  be  quite  interesting 
to  determine  if  this  adaptivity  could  be  employed  in  determining  the 
subregion  boundaries.  If  this  is  possible  then  variable  knots  splines 
could  become  a useful  image  segmentation  tool. 
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APPENDIX  A 


SOME  PROPERTIES  OF  NORMALIZED  B -SPLINES 


In  this  appendix  the  properties  of  the  normalized  B-splines  of 

order  k are  discussed.  The  discussion  will  be  limited  to  the  one 

dimensional  case  as  the  extension  to  a direct  product  of  splines  for 

two  dimensional  approximations  is  immediate. 

For  the  one  dimensional  case  a spline  S,  (x)  of  order  k with 

k,  N 

x 

N knots  approximating  f(x)  is  given  by 

Nx 

S,  ..  (x)  k 53  s.N.  ,r;x)  = f(x)  xe  fo , 1 1 
k,  N . . l l,  k — L J 

x 1 = 1 

where  N.  , (?,x)  are  the  normalized  B-splines  of  order  k satisfying 

t,  k 

the  following  recursion  relationship  [A-l]  over  the  knot  vector  % 
where 
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N +k 
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N.  , (ff;x)  = 
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• ,i  i-F-  1*k-1  ~ i + l . k-1 

i+k-1  l i+k  i + l 


(A-l) 


ot  herwi  se 


i’  i + k 


140 


TtI  s 


where  [T 1 is  as  follows 


[T]  = 


0 0 


3 0 


-6  3 


3 -3 


With  these  properties  we  can  now  demonstrate  an  interesting 
property  of  k*  ^ order  least  squares  splines.  A k order  least 

squares  spline  wit h N knot s for  f(x. ) i = 1 , 2 N (N^N  ) is  given 

by  the  solution  to  the  following  normal  equations: 


— N xl 

x 


■ tCNf  ff’tfxN 


(A-2) 


x y x 


where  f = [f(x  ) f(x  )1 

— Nxl  1 N 


If  we  take  N = k and  the  knot  vector  - to  consist  of  two  k order 
x 

knots  at  0 and  1 then  this  is  equivalent  to  determining  a k-1  degree 

least  squares  polynomial  for  fix.)  i = 1 N.  Alternately  if  we 

take  the  knot  vector  - to  contain  the  two  k order  knots  at  0 and  1 

and  a sufficient  number  of  internal  knots  so  that  N = N then  the 

x 

matrix  will  be  square  and  nonsingular  so  that  eq.  (A-2t 


becomes 
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1 xe 


Mi+1 


N ,(?;x)  = 
* » 1 


0 otherwise 


X)  = 1 


Note  that  eq.  (A-l)  indicates  that  Nx+k  knots  are  required  to  generate 
Nx  norrnalized  B-splines  of  order  k,  and  that  N , (F;x)  is  nonzero 
only  over  the  interval  [?.,  ?.+k>.  Also  a knot  may  have  multiplicity  p, 
up  to  k in  which  the  multiplicity  indicates  a discontinuity  in  the 
(k-(p+l ))  derivative  at  that  knot.  If  we  follow  Rice  [A-2]  and  adopt 
the  convention  that  the  spline  is  differentiable  of  order  0 or  -1  at  the 
knot  if  the  spline  is  continuous  or  has  a simple  jump  at  -q 
respectively.  a fourth  order  order  knot  at  F for  a cubic  or  4*  h 
order  spline  indicates  a simple  jump  in  the  spline  at  Figures  A.  1 
and  A.  2 illustrate  the  normalized  B-splines  of  order  4 for  the  knot 


vect  ors 


ll  = [0.  0.0,0,.  25,.  5,.  75,  1.  1,  1,  if 


12  = [0.  0,0,0,  .7.  . 8,  .9,  1.  1,  1,  ll1 

respectively.  corresponds  to  the  uniform  knot  case  while  * 
illustrates  the  effect  of  shifting  the  internal  knots  towards  1. 

I able  A.  1 lists  the  knots  over  which  N.  (*;x)  is  non-zero  for  i = 


1,2 


7.  Note  also  that  there  are  a total  of  11  knots  to  define  these 


7 nonzero  normalized  B-splines. 

Table  A.  1 and  figure  A.  1 serve  to  illustrate  the  relationship 

between  multiple  knots  and  the  differentiability  of  the  normalized 

B-splines.  Note  that  N^  ;x)  involves  a fourth  order  knot  at  x =0 

so  that  k-(p+l ) is  equal  to  -1  and  Nj  »x'  possesses  a simple 

jump  at  x = 0.  The  third  order  knot  at  x = 0 for  N1  A-,  ;x)  results  in 

2,4—1 

k-(p+l)  equaling  0 and  from  figure  A.  1 it  is  clear  that  N_  ("  ;x)  is 

l.  , 4 — “ 1 

merely  continuous  at  x = 0,  The  second  and  first  order  knots  at  x = 0 

for  N ("  ;x)  and  N.  A* . ;x)  respectively  result  in  N (- ;x)  being 
.5,41  4 , 4 —l  4 “1 

once  continuously  differentiable  and  N^  being  twice  continuous- 

ly differentiable  at  x=0.  The  same  sequence  of  events  is  true  for 

N5  4(Il;x)>  Nc,  4— l’X'  anfl  N7  4^— l’x'  at  x = 1 • 

Another  interesting  property  of  the  k order  normalized 

B-splines  arises  in  the  case  when  the  knot  vector  consists  of  solely 

two  k ^ order  knots  at  x = 0 and  x = 1.  In  this  case  the  k ^ order 


spline  with  N^  = k knots  can  be  shown  to  be  equivalent  to  a k-1  degree 

polynomial.  For  k = 4 the  four  nonzero  normalized  B -splines  are 

listedinTableA.2  fortheknot  vector  =(0,0,  0,0,  1,1,  1,1),  In 

4 

this  case  if  S,  ,(x)  = / S.N.  ,(c';x)  is  the  fourth  order  spline  with 
4,  4 ^ i i,4  — 

four  knots  and  ^|a. x' ~ ' is  a cubic  polynomial  it  is  simple  to  show 
that  a = (a^.a^.a^.a^)1  and  S = (Sj , S^,  S^,  S^)  are  related  through 
the  matrix  [T"|  by 


S = [N  f'‘lL 

corresponding  to  interpolation.  Thus  by  varying  the  number  of 

internal  knots  from  0 to  N-k  so  that  N = N k+k  = N,  we  can  deter- 

x 

mine  a least  squares  spline  that  goes  from  a least  squares  polynomial 
of  degree  k-1  to  an  interpolating  spline  of  order  k. 

In  summary  this  appendix  has  been  presented  to  illustrate  some 
of  the  computational  properties  of  the  normalized  R-splines.  This  has 
been  done  in  the  hope  that  the  interested  reader  might  be  better  able 
to  appreciate  them  from  a computational  standpoint  and  to  more 
effectively  apply  them  to  approximation  problems.  It  is  also  felt 
that  these  facts  coupled  with  the  recursion  formula  of  eq.  (A-l) 
should  give  the  reader  enough  information  to  implement  his  own 
spline  subroutines  with  variable  knots  without  an  excessive  amount 
of  difficult  y. 


14-1 


■Mi 


Figure  A.  1 (continued)  Normalized  4ih  Order 
B- Splines  For  Knot  Vector 

A (0,  0,  0,  0,  . 25,  . 5,  . 75,  1,  1,  1,  1) 


0 .1.2.3  .4.5.6. 
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